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AN  INVESTIGATION*  OF  THE  CHARACTERISTIC  AND  INSTABILITY  BEHAVIOR  OF  AN  AXIAL 
COMPRESSOR  IN  A  TURBOJET  ENGINE 


Inlet  Distortion  Research  Group,  Northwest  Polytechnical  University 


In  this  paper,  experimental  results  and  analysis  about  the  effect  of  the  inlet  pres¬ 
sure  distortion,  decrease  in  the  first  stage  turbine  nozzle  area  and  engine  test  run  program 
on  the  characteristic  and  instability  behavior  of  an  axial  compressor  in  a  turbojet  engine 
are  presented.  The  experiments  have  been  conducted  on  the  test  bed  on  a  turbojet  engine 
•with  a  9-stage  compressor,  whose  first  stage  is  transonic. 

The  results  show  that:  (1)  inlet  distortion  and  decrease  in  turbine  nozzle  area  may 
shift  considerably  the  speed  lines  of  compressor  leftward  and  cause  it  to  become  more  flat 
in  shape;  (2)  the  different  run  program  for  engine  test  may  produce  different  speed  line 
of  compressor  both  in  position  and  shape;  (3)  at  low  and  medium  speed,  the  instability  of 
compressor  exhibits  itself  as  a  rotating  stall,  and  at  high  speed,  a  surge,  whose  mode  de¬ 
pends  on  the  degree  of  inlet  distortion  and  decrease  in  turbine  nozzle  area. 


I .  FOREWORD 

The  characteristic  and  instability  behavior  of  a  multistage  axial 
compressor  is  a  problem  of  concern  to  design  personnel,  who  not  only  have 
to  understand  the  characteristics  of  the  axial  compressor  components  during 
the  tests,  but  more  important  is  to  understand  the  characteristics  of  the 
propulsion  installation  of  the  axial  compressor  or  characteristics  of  the  gas 

♦This  paper  was  read  at  the  Third  All  China  Engineering  Thermophysics  Conference 
at  Guilin  in  April  1980.  1 


turbine  installation  in  order  to  predict  and  analyze  the  matching  among 
components  and  characteristics  of  the  entire  machine.  Besides,  the  design 
personnel  require  to  understand  the  instability  behavior  of  a  compressor 
at  the  stall  line  in  order  to  extensively  study  the  cause  of  instability 
of  an  axial  compressor  at  the  stall  line  in  order  to  extensively  study 
the  cause  of  instability  of  an  axial  compressor,  prediction  of  unsteady 
vibration,  and  prevention  of  unsteady  vibration.  This  test  was  undergone 
on  a  turbojet  engine.  The  paper  stresses  the  introduction  of  influence 
of  various  factors  on  characteristics  of  an  engine  compressor  and  on 
instability  behavior.  The  testing  revolution  speed  was  run  between  87  and 
100  percent  of  the  rated  revolution  speed.  There  are  two  types  (of  anomaly 
meshes):  one  layer  of  36  eyes,  and  three  layers  of  36  eyes.  There  are  two 
types  of  first  stage  nozzle  area:  one  is  the  normal  area  at  the  design 
stage,  and  the  other  is  reduced  by  12  percent.  In  addition,  according  to 
different  starting  sequences,  curves  of  equal  revolution  speed  of  88.4  percent 
of  the  designed  revolution  speed,  and  the  operating  points  at  87,  90  and 
93.3  percent  of  the  designed  speed  were  measured. 

II.  EXPERIMENTAL  INSTALLATION  AND  EXPERIMENTAL  METHOD 

Figure  1  shows  the  test  installation  and  schematic  diagram  for  measure¬ 
ment  points.  In  front  of  the  engine  inlet  at  a  distance  of  about  0.5  of  the 
diameter  length,  a  rotatable  metal  mesh  is  installed  to  produce  the  required 
square  wave  anomaly  spectrum.  A  movable  tail  cone  was  installed  behind 
the  nozzles  of  the  engine  in  order  to  induce  unsteady  engine  vibration. 

While  studying  the  effect  of  inlet  anomaly  and  variation  of  the  first 
stage  nozzle  area  of  the  turbine,  the  engine  revolution  speed  was  maintained 
constant.  The  nozzle  openings  were  changed  and  the  tail  cone  moved  to  let 
the  operating  points  of  the  engine  move  at  a  constant  revolution  speed  curve 
until  unsteady  vibration  or  stall  appeared. 

In  studying  the  influence  of  the  starting  sequence,  starting  was 
conducted  according  to  the  following  two  sequences:  one  is  accelerating  from 
idling  revolution  speed  to  testing  revolution  speed  before  throttling  along 
the  equal  revolution  speed  curve.  Another  is  decelerating  from  the  designed 
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revolution  speed  (the  nozzle  opening  is  in  the  accelerating  state)  to  the 
testing  revolution  speed  before  throttling  along  the  equal  revolution  speed 
curve. 

III.  RESULTS  AND  ANALYSIS 

(A)  Influence  due  to  anomaly  of  inlet  gas  pressure 


Fig.  1.  Schematic  diagram  of  testing  installation 
and  testing  point  arrangement. 

Key:  (a)  Rotatable  mesh;  (b)  Movable  tail  cone; 

(c)  Total  pressure  in  steady  state;  (d)  Target  of 
total  pressure;  (e)  Total  pressure  in  dynamic  state; 

(f)  Overall  temperature;  (g)  Static  pressure;  (h) 

Dynamic  state;  (i)  Nine  point  target  of  total  pressure. 

While  Figs.  2,  3  and  4  are,  respectively,  n=0.9  and  one  layer  with  36 
eye  mesh  (covering  circumferential  angle  at  180°),  the  pressure  spectrum, 
total  pressure  circumferential  spectrum  and  wall  surface  static  pressure 
spectrum  at  the  inlet  cross  section  of  the  engine,  the  damping  function  of 
the  mesh  creates  the  circumferential  anomaly  spectrum  of  total  pressure 
nearing  a  square  wave;  however,  the  suction  function  on  a  high  and  low  pres¬ 
sure  region  gas  stream  by  the  compressor  causes  the  static  pressure  anomaly 
as  shown  in  Fig.  4. 

Figure  5  shows  the  influence  of  engine  compressor  characteristics  on 
the  inlet  gas  anomaly  while  n«0.9.  In  the  diagram,  three  equal  revolution 
curves  are  shown  for  cases  without  mesh,  one  layer  with  36  eye  mesh  (180° 
covering),  and  three  layer  mesh  with  36  eyes  (180®  covering).  The  dot  and 
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dash  curve  shows  the  results  of  the  compressor  component  test  on  a  platform. 
In  the  diagram,  three  curves  of  equal  nozzle  area  show,  respectively,  nozzle 
opening  positions  at  acceleration  (4>),  rated  power  (H) ,  and  the  maximum 
overload  (M)  condition. 


Fig.  2.  Pressure  spectrum  at 
engine  inlet. 


Fig.  3.  Total  pressure  circumferential 
spectrum  at  engine  inlet. 

Key:  (a)  Circumferential  angle;  (b) 

Anomaly  mesh. 


Fig.  4.  Static  pressure  circumfer-  Fig.  5.  Influence  on  engine  compressor 
ential  spectrum  at  engine  inlet.  characteristics  of  inlet  gas  anomaly. 

Key:  (a)  Circumferential  angle;  Key:  (a)  On  engine;  (b)  Component  test; 

(b)  Anomaly  mesh.  (c)  Point  of  unsteady  vibration;  (d)  Measure¬ 

ment  point;  (e)  Stall  line  of  component; 

(f)  Without  mesh. 


4 


As  revealed  by  the  test  results,  the  inlet  gas  anomaly  moves  leftward 
and  levels  the  equal  revolution  speed  curve;  the  greater  the  anomaly,  the 
more  pronounced  is  the  variation.  The  main  cause  of  "leftward  motion"  and 
"leveling"  is  possibly  weakening  of  the  toal  circulating  capacity  of  the 
compressor  with  reduction  in  total  flow  and  intensifying  of  perturbation 
between  stages  (lowering  of  efficiency  and  reduction  of  pressure  ratio). 

The  experimental  results  also  reveal  that  the  inlet  gas  anomaly  reduces 
the  allowance  of  unsteady  vibration  as  shown  in  Fig.  5.  Owing  to  leftward 
motion  and  leveling  of  the  equal  revolution  speed  curve,  the  lowering 
allowance  of  unsteady  vibration  is  the  overall  result  of  simultaneous  lowering 
of  pressure  ratio  and  flow. 

Figure  6  gives  the  characteristics  of  two  (anomaly  and  "net")  sub¬ 
compressors  and  the  average  characteristics  of  the  whole  compressor  set.  It 
is  apparent  from  the  figure  that  before  unsteady  vibration  of  the  compressor, 
the  pressure  ratio  of  the  anomaly  subcompressor  exceeds  the  pressure  ratio 
of  stability  state  unsteady  vibration.  This  is  because  in  an  anomaly  flow 
field,  the  rotor  blades  are  undergoing  periodic  variation  of  the  gas  stream 
attack  angle.  Under  this  instability  state,  the  proved  phenomenon  of  lagging 
dynamic  stall  (larger  attack  angle  can  be  borne  by  the  blades  than  in  the 
stability  state  without  stall)  to  let  the  operating  pressure  ratio  of  the 
anomaly  subcompressor  exceed  the  stability  state  pressure  ratio  (of  unsteady 
vibration)  while  the  compressor  still  maintains  a  steady  operation. 

Figures  7  and  8  show,  respectively,  the  situation  of  engine  unsteady 
vibration  at  two  different  anomalies  for  the  time  duration  of  compressor  exit 
total  pressure  p2*  and  the  static  pressure  ps  of  the  flow  measurement  sector. 
It  is  apparent  by  comparing  two  diagrams  that  the  anomaly  degree  has  an  effect 
on  the  unsteady  state  of  the  gas  stream.  It  can  be  considered  that  the 
"classical  unsteady  vibration"  appears  for  three  layer  36  eye  mesh,  and 
"deep  unsteady  vibration"  appears  for  one  layer  36  eye  mesh.  Figure  9  (with 
operating  condition  the  same  as  that  of  Fig.  7)  shows  the  photograph  of  a 
mesh  having  fallen  down  by  reverse  shock  of  the  shock  wave  during  unsteady 
vibration  of  the  engine. 


Fig.  6.  Compressor  characteristics  Fig.  7.  Time  duration  of  p£*  and  p  with 
of  anomaly  engine  unsteady  vibration  of  the  engine  atglx36  eye 

Key:  (a)  Point  of  unsteady  vibra-  mesh. 

tion;  (b)  Measurement  point;  (c)  Key:  (a)  Second;  (b)  Time. 

Component;  (d)  Anomaly  subcompres¬ 
sor;  (e)  "Net"  subcompressor; 

(f)  Average. 


Fig.  8.  Time  durations  of  p,*,  p  Fig-  9-  Photograph  of  a  mesh  falling  down, 
and  T  *  for  unsteady  vibration  of8  by  a  shock  wave  during  unsteady  vibration, 
the  engine  with  a  3x36  eye  mesh. 

Key:  (a)  Time;  (b)  Second. 
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(B)  Influence  of  first  stage  nozzle  area  of  the  turbine. 

Figure  10  shows  the  variation  of  engine  compressor  characteristics  when 
the  first  stage  nozzle  area  of  a  turbine  is  reduced  by  12  percent.  Reduction 
of  first  stage  nozzle  area  moves  the  operating  points  of  the  engine  upward; 
the  compressor  equal  revolution  speed  lines  move  leftward  and  become  level; 
the  lower  the  revolution  speed,  the  more  level  are  the  equal  revolution  speed 
curves,  even  appearing  with  the  left  branch  of  the  revolution  speed  curves 
at  slightly  positive  values  of  the  inclination.  Figure  11  shows  the  time 
durations  of  P2*>  P  and  T_*  (with  n=0.87)  undergoing  unsteady  vibration  when 
the  first  stage  nozzle  area  is  reduced  by  12  percent.  It  is  apparent  from 
the  diagram  that  the  shock  wave  produced  during  unsteady  vibration  arrives 
at  the  flow  measurement  sector  within  approximately  6  to  8  milliseconds.  At 
that  time,  a  peak  value  of  p^  appears.  According  to  analysis,  the  unsteady 
state  is  the  "deep  unsteady  vibration". 

(C)  Influence  of  starting  sequence 

Figure  12  shows  the  equal  revolution  speed  curves  of  n=0.884  measured 
at  different  starting  sequence  and  the  operating  points  of  n=0.933,  0.9  and 
0.87.  The  small  black  circles  in  the  diagram  show  the  measurement  values 
from  idling  revolution  speed  accelerating  to  the  testing  revolution  speed; 
the  small  white  circles  represent  the  measurement  values  from  design  of 
revolution  speed  expanding  nozzle  area  decelerating  to  the  testing  revolu¬ 
tion  speed.  From  the  limited  data  shown  in  the  diagram,  it  is  apparent  that 
we  can  see  the  "multiple"  characteristics  of  the  equal  revolution  speed  curves 
of  the  axial  compressor;  i.e.,  different  equal  revolution  speed  curves  can 
be  obtained  with  different  manners  of  starting .  The  appearance  of  "multiple" 
characteristics  may  possibly  be  related  to  the  "lagging  phenomenon"  of 
rotational  stall  and  interstage  perturbation  of  the  axial  compressor. 

(D)  Influence  of  engine  revolution  speed 

The  engine  revolution  speed  is  one  of  the  major  factors  affecting  the 
unsteady  state  of  a  compressor.  As  revealed  by  considerable  analyses  and 
experimental  studies,  during  low  speed  generally  rotational  stall  appears  on 
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the  engine  at  the  stabilized  boundary.  However,  during  high  revolution  speed, 
unsteady  vibration  frequently  occurs.  This  study  also  proves  this  point. 


Fig.  10.  Influence  of  compressor  charac¬ 
teristics  by  reducing  12  percent  of  the 
first  stage  nozzle  area  while  n=0.87. 

Key:  (a)  On  engine;  (b)  Component;  (c) 

Small  nozzles  for  the  first  stage;  (d) 
Normal  nozzles  for  the  first  stage;  (e) 
Point  of  unsteady  vibration;  (f)  Stall 
curve. 


•  JiiiiZSri.i 


Fig.  11.  Time  durations  of  p-*, 
p^  and  T2*  during  unsteady  vibra¬ 
tion  while  the  first  stage  nozzle 
area  is  reduced  by  12  percent. 
Key:  (a)  Time;  (b)  Second. 


Fig.  12.  Influence  on  engine  compressor 
characteristics  of  starting  sequence. 

Key:  (a)  Component  test;  (b)  Decelerating 
process;  (c)  Accelerating  process;  (d) 
g Operating  curve. 


Figure  13  shows  the  behavior  overtime  of  o'j,  T£  and  pg 
and  interstage  p.^  j^and  p*  3  when  n  =  0.69  (called  "upper 
breathing  point)  at  stall;  the  frequency  is  46  percent  of 
the  rotational  speed.  It  is  so-called  "lower  point  of  un¬ 
steady  vibration")'.  However,  when  n=0.87,  the  unsteady  state 
appeared  as  "deep  unsteady  vibration" shown  in  Fig.  11. 


Figure  13.  The  behavior  overtime 
of  pg,  PL>1/  P]_.3'  P2'  and  T 2  When 
n=0.69  and  at  stall. 

Key;  (*)  Time. 

The  problem  of  unsteady  vibration  of  a  compressor  can  be  analyzed  by 

using  the  Helmholtz  acoustical  resonator  model.  During  unsteady  vibration 

of  the  compressor,  the  gas  stream  vibrates  at  the  front  and  rear;  while  using 

a  resonator,  the  gas  column  vibration  at  the  neck  portion  becomes  a  model. 

The  motion  equation  of  the  gas  column  at  the  neck  portion  of  the  resonator 

can  be  written  as  [31:  j?x  /  j  \ 

J  ~  +  at[-£s-\x  -  0 

4*  \  V,L') 

Obviously,  the  vibration  frequency  (of  the  gas  column) 

«  "■  AJVtL r. 

In  the  equation,  a  is  the  sound  speed;  Ac  is  the  equivalent  area  of  the 
compressor  ring  surface;  is  the  volume  behind  the  compressor;  and  Lc 
is  the  equivalent  length  of  the  compressor. 
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B  -  V/1*L ,  -  (l//2«)[vr  Vt/A,L.  1 

Through  experiments  and  analyses,  Greit2er  [4,  5]  proved  the  existence  of  a 
critical  value  Bcritical*  When  B<Bcritical*  rotat*onal  stall  appears  at  the 
steady  boundary;  unsteady  vibration  appears  when  B  ^  Bcriticai •  I*  is  apparent 
from  the  above  equation  that  the  value  of  B  is  determined  by  the  rotational 
speed  for  a  machine  of  a  definite  structure.  Therefore,  unsteady  vibration 
easily  occurs  at  high  revolution  speed.  Rotational  stall  often  appears 
during  low  rotational  speed.  In  the  result  of  this  experiment,  it  is  the 
abrupt  type  stall  while  n^  0.69  (B^  0.79);  it  is  the  unsteady  vibration  when 
n>  0.87  (B^  0.98). 


It  should  be  pointed  out  that  by  using  the  Helmholtz  resonator  model, 
the  actual  system  of  the  influences  on  the  unsteady  state  by  the  following 
factors  is  not  considered:  position  of  throttling  valve,  inlet  gas  anomaly, 
and  operating  point  positions  of  compressor.  However,  for  rotational  speed, 
the  revelation  trend  of  the  above  mentioned  expression  equations  is  consistent 
with  the  experimental  results. 


IV.  CONCLUSIONS 


(1)  For  an  axial  compressor  inside  a  turbojet  engine,  the  characteristics 
are  not  only  determined  by  the  characteristics  of  the  compressor  components, 
but  are  also  related  to  the  structure  layout  and  operating  conditions  of 
other  components  in  the  compressor  system.  The  inlet  gas  anomaly  and  reduction 
of  the  first  stage  nozzle  area  of  the  turbojet  engine  can  move  leftward  and 
become  more  level  for  the  equal  revolution  speed  curves;  the  lower  the 
revolution  speed,  the  more  pronounced  is  the  leftward  motion  and  leveling. 

A  different  starting  sequence  also  will  cause  variation  to  the  shape  and  posi¬ 
tion  of  the  equal  revolution  speed  curves. 

(2)  The  unsteady  state  appearing  at  the  steady  boundary  of  the  turbojet 
engine  compressor  is  mainly  determined  by  volume  and  revolution  speed  of  the 
combustion  chamber.  At  the  same  time,  the  unsteady  state  is  related  to 
operating  conditions  and  the  system  of  the  compressor.  At  n^0.69  for  the 
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turbojet  engine  in  the  experiment,  the  appearing  unsteady  state  is  the  abrupt 
type  stall;  during  n  ^0.87,  it  is  the  unsteady  vibration.  In  these  situations, 
"classical  unsteady  vibration"  appears  when  three  layers  of  36  eye  mesh 
appear.  While  the  first  stage  nozzle  area  of  the  turbojet  engine  is  reduced 
by  12  percent,  the  "deep  unsteady  vibration"  appears. 

(3)  Due  to  lagging  of  the  dynamic  stall,  the  operating  pressure  ratio 
of  the  subcompressor  during  low  inlet  gas  pressure  may  exceed  the  pressure 
ratio  of  stability  state  unsteady  vibration  of  a  homogeneous  inlet  gas. 

The  experiment  of  this  paper  was  collectively  completed  by  the  Inlet 
Distortion  Research  Group. 
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PREDICTION*  OF  ONSET  OF  ROTATING  STALL  USING  SMALL  PERTURBATION  THEORY  FOR 
AXIAL  FLOW  COMPRESSORS 

Liu  Zhiwei** 

Northwest  Polytechnical  University 

This  paper  analyses,  generalizes  and  compares  different  onset  conditions  of  the 
predicting  rotating  stall,  which  are  based  on  small  perturbation  theory.  They  are 
correlated  to  experimental  data  of  rotors  for  single  stage  axial  flow  compressors.  Good 
agreement  between  theoretical  values  and  experimental  values  is  found.  Mechanism 
for  the  onset  condition  of  rotating  stall  is  analysed.  Effectiveness  of  methods  of 
prediction  is  proved.  The  extent  of  application  of  methods  is  extended. 

I .  FOREWORD 

The  prediction  of  onset  of  rotating  stall  is  related  to  the  determination 
of  the  allowance  index  and  compressor  stall  boundary.  Up  to  now,  there  is 
still  no  method  to  be  followed.  Starting  from  different  viewpoints,  quite 
a  few  sources  of  information  proposed  the  basis  of  evaluating  the  methods. 

It  seems  that  those  conclusions  do  not  have  generalized  significance  owing 
to  lacking  of  more  experimental  data.  Or,  it  is  difficult  to  more  accurately 

*The  paper  was  read  at  the  Third  All  China  Engineering  Thermophysics  Confer¬ 
ence  at  Guilin  in  April  1980. 

**Also  taking  part  in  the  experimental  studies  of  the  paper  were  Zhang  Zhang- 
sheng,  Zhang  Weide,  Shi  Jingxun,  Huang  Jiangong,  Wang  Zongyuan,  and  Lin  Qixun. 
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determine  the  empirical  formula  because  of  too  great  a  range  of  the  empirical 
data.  It  seems  that  the  criterion  established  by  analyzing  the  stability  of 
fluid  passing  through  the  blade  array  on  the  basis  of  small  perturbation 
theory  can  basically  get  to  the  nature  of  the  problem.  In  addition,  good 
agreement  between  the  theoretical  and  experimental  values  can  be  obtained. 

The  work  of  this  aspect  has  been  proven  by  G.  R.  Ludwig  et  al  [1]. 

Noticing  different  premises  and  treatment  methods  of  several  small 
perturbation  theories  in  analyzing  the  fluid  stability,  there  are  different 
criteria  in  predicting  stall  onset  conditions.  In  addition,  on  viewing 
the  data  in  [1]  correlating  only  at  a  low  M  number  (smaller  than  0.2)  and 
large  hub  ratios,  the  influences  of  fluid  compressibility  and  three-dimensional 
effects  are  considerably  different  than  for  rotors  used  in  practice.  In  order 
to  discuss  the  effectiveness  in  practical  applications  of  the  above-mentioned 
theoretically  obtained  criteria,  in  this  paper  correlation  is  made  between 
rotor  data  and  theory  for  comparing  various  types  of  small  perturbation 
theories  in  analyzing  the  stall  onset  points.  This  goes  a  further  step  in 
proving  the  effectiveness  of  predicting  the  onset  of  a  rotating  stall  of  a 
single  blade  array;  the  range  of  application  is  extended. 

II.  ANALYSIS  OF  SMALL  PERTURBATION  THEORY 

It  is  assumed  that  the  passing  through  the  blade  array  by  the  gas  stream 
is  two-dimensional,  incompressible  and  unsteady,  and  the  blade  array  is 
substituted  by  a  model  of  an  excitation  disc  as  shown  in  Fig.  1.  At  the 
same  time,  it  is  assumed  that  the  perturbation  amount  of  the  gas  stream  can 
be  neglected  compared  to  the  corresponding  average  gas  flow.  This  is  the 
assumption  of  small  perturbation. 

From  the  method  of  solving  for  the  perturbation  amount  of  upper  and  lower 
streams  of  a  flow  field  of  the  blade  array,  the  small  perturbation  theory  can 
be  divided  into  two  categories:  one  is  the  series  form  solution  in  perturbation 
speed  potential  function  at  the  upper  stream  of  the  blade  array  in  order  to 
obtain  the  onset  conditions  of  evaluating  the  rotating  stall  [2).  Another 
category  is  the  multiple  index  form  flow  function  at  the  upper  and  lower 
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stream  blade  array  to  obtain  the  onset  conditions  of  evaluating  the  rotating 
stall  [1,  3,  4]. 


In  order  to  obtain  the  solution  of  the  second  category  flow  function, 
in  this  paper  a  motion  equation  of  perturbation  vorticity  describes  the  flow 
fields  at  the  upper  and  lower  streams  of  the  blade  array. 

&  (11  -I.  -5SL  j.  (V  J-  _  n  /i  \ 


+  a/  +  *)!SL  +  0'  +  _o  (i) 
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Under  the  given  onset  conditions  and  boundary  conditions,  the  perturbation 
vorticity  expressed  in  the  form  of  a  multiple  index  can  be  obtained: 


C  —  Co* 


By  relating  the  perturbation  flow  function  and  perturbation  vorticity 
satisfying  the  continuity  equation,  we  get: 

W  -  -  f  (3) 


Substitute  eq.  (2)  into  eq.  (3)  and  after  variable  replacement,  the 
solution  in  the  form  of  Poisson's  equation  can  be  derived,  in  obtaining  the 
expression  equation  with  the  same  form  as  the  perturbation  flow  function  in 
[3]: 

“  Ai*  +  Bjt 

U.  (-S-KIJJ 


If  the  same  coordinating  and  boundary  conditions  are  used  as  in  [1],  the 
homogeneous  equation  set  of  the  unknown  constants  A^,  and  can  be  obtained. 
According  to  conditions  of  non-singular  solution  of  this  equation  set,  one 
root  of  the  physical  significance  of  the  characteristic  equation  can  be 
determined,  thus  obtaining  the  onset  conditions  of  a  rotating  stall  as  pointed 
out  in  [1] : 

When  p2(®»  y) -constant,  the  damped  factor  is: 
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And  when  p2(0  ,  y) “constant,  the  damped  factor  is: 

TSi'-JL-\\  +  X  +  Si  -  S'&S,  -  ^  (l  +S?)X'l  (6) 

nil  2  +  f  l  2  * 

If  Cj  >  0,  the  gas  stream  is  stabilized;  if  <0,  the  gas  stream  is  not 
stabilized;  if  Cj*0,  this  is  the  neutral  stabilization  condition.  From  eq.  (5) 
or  eq.  (6),  the  onset  point  of  a  rotating  stall  can  be  determined. 

For  comparison,  conclusions  of  various  methods  should  be  homogeneous. 

This  paper  conducts  further  analysis  of  the  onset  conditions  in  [3]  and  [4]. 


From  [3],  the  onset  condition  is: 

f'  -  2(|  +  crf^/e if 

In  addition,  it  is  noted  that  £  »  X(1  + 
1  +  X  +  Si  -  [S,(l  +  5?)*V2]  -  0 


and  r-0f/actgA. 


If  the  influence  of  time  lag  of  the  adhering  layer  is  neglected,  and 
from  [4]  it  is  determined  that  the  increment  rate  is  equal  to  zero  to  determine 
the  stall  onset  point,  similarly,  the  same  expression  equation  as  eq.  (8) 
can  be  obtained,  thus  determining  that  there  is  no  inherent  difference  in 
the  linear  theoretical  part. 

Comparing  eq.  (5)  or  eq.  (6)  with  eq.  (8),  it  is  apparent  that  eq.  (8) 
is  the  special  case  of  the  two  other  equations.  If  the  variation  of  exhaust 
angle  with  gas  inlet  angle  is  not  considered,  these  equations  have  the  same 
form.  The  basic  cause  of  their  differences  is  the  diversion  characteristics 
of  the  blade  array. 

III.  SOLVING  EXPERIMENTAL  PROBLEMS 


i 
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Using  correlated  comparison  between  the  analysis  results  and  four  rotor 
data,  there  are  the  isolated  rotor  data  (variable  energy  design  of  the 
rotor  along  the  blade  height)  and  the  rotor  data  provided  in  [5-7]  of  a  medium 
hub  ratio  of  the  simultaneous  experiment.  In  correlation,  influence  on  the 
rotor  is  not  considered.  The  related  parameters  of  these  test  pieces  are 
compiled  in  Table  1. 

Table  1.  Rotor  parameters 


- - ~~ _  No. 

•  — 

SXfll 

[5] 

[7] 

(cjE  tfc 

■eh 

1.20 

1.20 

1.53 

****«««»  (d) 

0.87 

0.88 

1.2 

(ep  «  tfc 

0.61 

0.46 

0.42 

0.40 

«  «  tfc(f) 

1.94 

1.78 

3.28 

6.28 

48.2® 

,  72.8® 

62.3* 

53.6® 

n  A  (h) 

42.6® 

39.5° 

22.3® 

21.6® 

(IP  * 

1.17 

0.998 

1.623 

1.881 

*  All  parameters  related  with  the  cross  section  mean  the  arithmetic 
mean  parameters  for  an  average  radius. 

Key:  (a)  Parameters;  (b)  Author's  experiments;  (c)  Pressure  ratio; 

(d)  Relative  M  number  of  incoming  flow  at  the  blade  tip;  (e)  Hub 
ratio;  (f)  Aspect  ratio;  (g)  Installation  angle*;  (h)  Bending  angle; 
(i)  Viscosity. 


Data  compilation  method:  The  average  parameters  of  a  gas  stream  passing 
through  the  rotors  are  substituted  by  mass  average  parameters.  The  geometrical 
characteristics  cross  section  of  rotors  is  selected  for  the  average  radius  of 
the  area.  The  data  compilation  progresses  on  the  assumption  of  cylinder  flow. 
The  loss  characteristic  of  the  rotor  blade  array  is  calculated  according  to 
the  following  equation. 

rDKri.  -  (?) 

The  diversion  characteristic  is  obtained  by  averaging  the  measured  distribution 
parameters  along  the  blade  height.  Figure  2  lists  the  average  characteristics 
of  subsonic  rotors  in  the  author's  experiments  and  the  transonic  rotor  provided 
by  [7].  There  is  no  repetition  for  characteristics  of  other  rotors  because 
of  similarity. 


Processing  of  experimental  data:  In  Fig.  2,  the  inlet  gas  angle 
is  reduced  to  the  flow  pattern  that  the  rotating  stall  will  immediately 
appear;  however,  flow  is  still  maintained  at  the  corresponding  point  of 
normal  stabilized  flow,  i.e.,  the  so-called  neutral  stabilized  boundary 
point.  This  is  the  stall  boundary  point  as  it  is  usually  meant.  The  inlet 
gas  angle  corresponding  to  this  point  is  the  inlet  gas  angle  of  the  onset 
point  of  a  rotating  stall  determined  by  the  experiment.  For  an  abrupt  type 
rotating  stall,  the  blade  array  characteristic  curve  after  the  neutral 
stabilized  boundary  point  will  have  interruption  point.  After  the  curve 
passes  through  the  interruption  point,  the  natural  extension  is  the  approximate 
processing  method  in  data  correlation.  Data  processing  is  done  according 
to  the  method  of  least  squares  for  the  characteristics  of  various  rotors.  In 
addition,  the  following  polynomial  is  taken  for  curve  simulation. 


G 
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Fig.  2.  Average  characteristics  of  subsonic 
and  transonic  rotors. 

Key:  (a)  Symbol;  (b)  Author's  experiments; 
(c)  According  to  [7]. 


IV.  RESULTS  OF  CORRELATION  BETWEEN  EXPERIMENTS  AND  THEORY 


The  rotor  characteristics  after  processing  are  substituted  into  eq.  (5) 
or  (6),  and  eq.  (8).  The  corresponding  inlet  gas  angle  at  onset  of  the 
theoretical  rotating  stall  is  obtained;  comparison  is  made  between  the 
theoretical  and  experimental  values.  Figure  3  lists  curves  of  the  damped 
factors  rCj/nU  (varying  with  Sj)  calculated  with  correlation  using  the  author's 
experiment  data  and  eq.  (5).  It  is  apparent  from  Fig.  3  that  good  agreement 
exists  between  the  theoretical  and  experimental  values. 


Fig.  3.  Theoretical  stability  characteristic 
[using  eq.  (5)]  of  rotors. 

Key:  (a)  Steady  region;  (b)  Unsteady  region; 

(c)  Onset  point  of  experiment. 


Tables  2  and  3  list,  respectively,  the  rotor  data  according  to  the  author's 
experiment  and  using  [5-7]  and  adopting  eq.  (5)  or  eq.  (6)  for  the  correlation 
calculation  in  obtaining  the  onset  inlet  gas  angle  and  the  corresponding 
experimental  value.  Table  4  lists  the  compared  results  of  two  rotors  with 
the  above-mentioned  (similar)  calculation  using  eq.  (8). 


From  comparison  of  four  rotors  listed  in  the  tables  mentioned  above,  it 
is  apparent  that:  Notwithstanding  the  considerable  difference  in  geometrical 
parameters  and  performance  characteristics  of  the  rotor,  there  is  good 
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Table  2.  Comparison  of  theoretical  value  [according  to  eq.  (5)  and 
eq.  (6)]  of  the  onset  point  inlet  gas  angle. 
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0.54 
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26.1 
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0.60 

1 

0.71 

27.3 

27.0 

27.2 

0.96 

0.35 

1 

0.84 

27.7 

20.9 

27.1 

2.9 

2.4 

1 

Key:  (a)  Experimental  value;  (b)  Theoretical  value;  (c)  Degree;  (d) 
Relative  error  percentage;  (e)  Number  of  stall  groups;  (f)  Constant; 
(g)  Experiment  and  theory. 


Table  3.  Comparison  of  theoretical  value  [according  to  eq.  (5)]  and 
experimental  value  of  the  onset  point  inlet  gas  angle. 
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0.8 
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34.6 

1.0 
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Key:  (a)  Degree;  (b)  Theoretical  value;  (c)  Experimental  value;  (d) 

According  to. 

Table  4.  Comparison  of  theoretical  value  [according  to  eq.  (8)]  of 
the  onset  point  inlet  gas  angle. 
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Key:  (a)  According  to  data  in  author's  experiments;  (b)  Degree;  (c) 
Theoretical  value;  (d)  Experimental  value;  Ce)  According  to. 
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agreement  between  the  inlet  gas  angle  at  stall  onset;  their  relative  error 
is  within  5  percent  except  in  individual  points  as  high  as  10  percent. 
Secondly,  differences  are  brought  about  by  the  onset  criterion  of  the  boundary 
condition  at  the  lower  stream  of  different  blade  array  outlets;  there  is  no 
major  influence  in  estimating  the  theoretical  values.  Thirdly,  as  revealed 
in  the  correlation  comparison  of  eq.  (5)  and  eq.  (8),  results  of  the  two  are 
satisfactory.  This  explains  that  the  diversion  characteristic  of  the  blade 
array  is  secondary  in  evaluation  of  the  onset  point;  however,  the  loss 
characteristics  have  a  determined  function;  this  has  a  restraining  effect 
on  the  diversion  characteristics.  In  other  words,  the  instability  factor 
of  the  blade  array  is  mainly  determined  by  the  related  term  of  the  loss 
variation  rate  X'.  The  term  related  to  the  exhaust  gas  angle  variation 
rate  82'  may  be  a  steady  factor,  and  may  be  an  unsteady  factor.  However, 
the  influence  on  stall  onset  point  is  very  slight,  even  to  a  negligible 
degree.  Attention  is  necessary  that  the  diversion  capability  of  the  blade 
array  and  the  loss  itself  created  are  the  factors  maintaining  the  stability. 

As  revealed  by  the  data  correlation  described  above,  although  there 
are  unavoidable  three-dimensional  influences  existing  among  various  rotor 
blade  arrays,  however,  this  three-dimensional  characteristic  does  not  play 
an  important  role  in  the  prediction  of  a  rotating  stall.  In  addition,  an 
average  blade  array  characteristic  can  be  found  as  representation  to  predict 
the  stall  onset  point.  It  is  also  indicated  that  the  fluid  compressibility 
has  a  secondary  position  in  the  influence  of  stall  onset  conditions.  This 
can  be  determined  by  going  a  further  step,  by  taking  into  account  the 
reliability  degree  and  adaptation  range  of  the  prediction  theory  of  two- 
dimensional  incompressible  unsteady  small  perturbation  of  fluid  inertia  in 
the  blade  grid  passage  but  neglecting  the  lagging  of  the  adhering  layer. 

It  should  be  point  out  that  the  theory  of  predicting  the  stall  onset 
mentioned  above  is  derived  from  the  appearance  of  a  rotating  stall  flow 
pattern.  However,  in  practice  often  not  a  single  unsteady  flow  pattern 
appears  along  the  blade  height  to  a  rotor  blade  array;  this  is  a  blended 
unsteady  flow  pattern  of  a  rotating  stall  and  unsteady  vibration.  In  this 
situation,  this  is  a  problem  worth  discussing  in  theoretically  predicting 
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how  effective  it  is.  On  this  point,  detailed  measurements  were  made  on 
dynamic  characteristics  of  unsteady  state  points  of  rotors  of  the  author's 
experiments.  The  results  are  shown  in  Figs.  4  and  5. 


Fig.  4.  Variation  of  circumferen¬ 
tial  time  averaged  total  pressure 
restoration  coefficient  along  the 
blade_height  at  the  rotor  inlet 
when  n=0.33. 

Key:  (a)  Symbol;  (b)  State;  (c) 

Before  stall;  (d)  After  stall; 

(e)  Throttling  again;  (f)  Complete 
unsteady  vibration. 


Fig.  5.  Dynamic  waveform  of  rotor  with 
appearance  of  unsteady  flow  pattern  when 
n=0.84. 

Key:  (a)  Pressure  transmitter;  (b) 

Second. 


Figure  4  qualitatively  explains  the  distribution  of  loss  of  inlet  energy 
along  the  blade  height  before  and  after  stall.  In  the  rotating  stall  region 
of  the  blade  height,  all  the  restoration  coefficients  of  the  total  pressure 
are  smaller  than  the  restoration  coefficients  of  total  pressure  during 
normal  flow.  However,  the  closer  to  the  portion  of  the  blade  tip  (r/r^  0.97), 
the  greater  the  restoration  coefficients.  This  explains  that  there  is 
energy  input  into  the  sector.  The  energy  comes  from  the  reverse  direction 
flow  of  the  outlet.  This  indicates  that  after  a  stall  there  exists  the 
blending  of  two  flow  patterns  (rotating  stall  and  unsteady  vibration)  along 
the  entire  blade  height.  At  the  same  time,  it  is  apparent  from  the  diagram 
that  with  increasing  throttling,  the  proportion  occupied  by  the  reverse 
direction  flow  region  along  the  blade  height  is  increased  until  the  complete 
disappearance  of  the  rotating  stall  and  the  flow  pattern  (of  unsteady  vibra¬ 
tion)  occupies  all  the  blade  height  region.  This  can  explain  (by  going  a 
further  step)  the  interrelated  situation  and  conversion  feature  of  the 
unsteady  flow  pattern.  21 


Figure  5  shows  that  the  low  frequency  oscillation  of  approximately  4  to  S 
hertz  is  that  of  an  unsteady  vibration  wave,  carrying  the  rotating  stall  wave 
with  frequency  approximately  at  120  hertz.  However,  the  first  appearance  of 
the  rotating  stall  wave  can  still  be  seen.  Approximately  0.4  second  later, 
a  low  frequency  unsteady  vibration  wave  then  appears.  As  revealed  by 
results  correlating  with  this  blending  unsteady  flow  pattern,  there  is  a 
good  agreement  between  the  theoretical  prediction  and  experimental  results. 

This  seems  to  explain  that  the  above  mentioned  linear  prediction  theory  can 
be  extended  to  the  evaluation  of  an  unsteady  flow  pattern. 

V.  CONCLUSIONS 

By  establishing  the  stability  criterion  of  analyzing  the  fluid  passing 
through  the  blade  array  on  the  basis  of  small  perturbations,  a  generalized 
expression  equation  (5)  or  (6)  can  be  derived.  Criteria  of  other  forms  can 
be  considered  as  special  cases.  By  analyzing  the  stall  onset  point  using 
the  two  dimensional,  incompressible  and  unsteady  flow  model,  it  can  be 
effectively  applied  to  the  correlation  of  actual  compressor  rotor  data; 
the  average  performance  can  be  used  as  the  basis  of  correlation.  The  three- 
dimensional  characteristic  and  compressibility  have  little  influence  on  the 
stall  onset  conditions.  In  the  regime  of  controlling  the  onset  of  a  rotating 
stall,  the  average  loss  characteristic  of  the  rotor  similarly  plays  a 
dominant  role.  To  the  appearing  situation  of  a  blending  unsteady  flow  pattern, 
the  above  mentioned  criterion  seems  feasible  to  be  applied.  More  experimental 
data  are  required  for  proof. 

As  revealed  by  the  above  mentioned  results,  it  is  possible  in  analysis 
of  the  criterion  basis  that  the  prediction  skill  of  compressor  unsteady  boundary 
and  stall  allowance  can  be  developed. 
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Key:  (a)  Polynomial  coefficient;  (b)  And;  (c)  Constants  in  solution  of  flow 

function;  (d)  Length  of  blade  chord;  (e)  Multiple  index  in  solutior  of  pertur¬ 
bation  vorticity;  (f)  Respectively  the  virtual  and  real  portions  of  C;  (g) 
Combination  parameter;  (h)  Unit  of  virtual  number;  (i)  Number  of  stall  groups; 
(j)  Percentage  of  design  revolution  speed;  (10  Static  pressure  of  gas  stream 
in  region  i;  (1)  Total  pressure  of  gas  stream  in  region  i;  (m)  Total  pressure 
of  gas  stream  relative  to  the  moving  blade  in  region  i;  (n)  Radius  of  an 
element  in  the  blade  array;  (0)  Cosecant  of  relative  inlet  gas  angle  of 
blade  array  in  region  i;  (p)  Cosecant  of  absolute  inlet  gas  angle  of  blade 
array  in  region  i;  (q)  Time  variable  relative  to  stationary  coordinates  of 
blades;  (r)  X  direction  component  of  perturbation  speed  of  gas  flow;  (s) 

X  direction  component  of  average  velocity  of  the  gas  stream;  (t)  Y  direction 
component  of  perturbation  speed  of  gas  stream;  (u)  Y  direction  component  of 
average  speed  of  gas  flow;  (v)  Loss  coefficient  of  total  pressure;  (w)  Coor¬ 
dinate  variables  fixing  on  blades:  x  along  the  axial  direction  and  y  along 
the  tangent  direction;  (x)  Relative  gas  flow  angle  (and  tangent  direction); 

(y)  Blade  installation  angle  (and  tangent  direction)  of  an  element  blade 
array;  (2)  Perturbation  vorticity  of  gas  stream;  (aa)  Loss  coefficient  of 
total  pressure  of  axial  component  definition  of  dynamic  pressure  head  at  inlet; 
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(Continuation  of  Key  of  a  table  in  the  previous  page] 

(ab)  Restoration  coefficient  of  total  pressure  at  inlet;  (ac)  Perturbation 
flow  coefficient  of  gas  stream;  (ad)  Subscripts;  (ae)  Region  of  gas  flow; 

(af)  Of  rotating  stall  onset  condition;  (ag)  Of  blade  tip;  (ah)  Inlet  region 
of  blade  array;  (ai)  Outlet  region  of  blade  array. 
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SOME  CONSIDERATIONS  CONCERNING  THE  PRINCIPAL  EQUATION  FOR  THE  AERODYNAMIC 
CALCULATION  OF  THE  NON-RADIAL  CALCULATION  STATION  ON  S_  STREAM-SURFACE  FLOW 
IN  TURBOMACHINES 


Zhang  Yujing 

Harbin  Marine  Boiler  and  Turbine  Research  Institute 
Received  10  October  1980 


Based  on  Wn’i  equations  expressed  in  terms  of  non-orthogonal  curvilinear  coor¬ 
dinates,  this  paper  derives  a  y-direction  equation  of  motion  for  the  aerodynamic 
calculation  of  the  non-radial  calculation  station  on  the  S«  stream-surface  flow  in  a 
turbomachine  blade  passage.  Thereby  it  is  proved  that  the  equation  applied  by 
Hetherington  (1974)  is  correct.  By  means  of  transformation,  it  is  pointed  out  that 
'he  principal  equation  (9)  of  this  paper  is  in  full  agreement  with  the  equation  (11) 

•■t  [2]. 


In  aerodynamic  calculation  of  the  S2  stream  surface  of  turbomachines, 
while  computation  stations  are  non-radial  direction  computation  stations  within 
the  gap  of  blade  arrays,  in  1974  Heterington  used  the  following  equation  (see 
equation  (30)  in  [1]): 

P  by  1  —  ML  l  WW*  2  by  r  l 

+  —  — —  [(\  —  U\)emt  +  A/’.tanftdac]  (I) 

r(l  — 


The  definition  of  the  (y,z)  coordinate  system  and  angles  e  and  6  used 
are  shown  in  Fig.  1*. 

Manuscript  [2]  expressed  disagreement  on  the  correctness  of  equation  (1). 


Fig.  1.  (y,z)  coordinate  system 

used  in  equation  (1). 


Starting  from  the  three-dimensional  flow  fundamental  equation  [3]  of 
turbomachines  with  arbitrary  curve  coordinates,  the  paper  derives  the  y 
direction  motion  equation  for  calculation  of  the  S2  flow  surface  of  the 
(y,  z,  <)>)  coordinate  system  shown  in  Fig.  1;  a  proof  was  given  that  equation 
(1)  is  correct.  By  transformation,  it  is  pointed  out  that  the  principal 
equation  (9)  derived  is  in  agreement  with  equation  (11)  of  manuscript  [2]. 

[ 

r 

! 

DERIVATION  OF  Y  DIRECTION  MOTION  EQUATION 


The  continuity  equation  and  motion  equation  of  an  S2  flow  surface  using 
an  arbitrary  curve  coordinate  is  [3] : 


(rpwW *u»  $u)  +  A  (tfiWW  <*u  sin0u)  m  0 

9xv  Ox1 

W 1  (  9m,  _  9mx\  Wj  9(Vtr)  _  9i_  _  _  9s_ 

dx>)  r'  9,'  9x'  9xx 


(2) 

(3) 


*The  direction  a  of  a  coordinate  in  the  paper  is  the  direction  z  of  manuscript 
[2].  Except  those  specified,  the  other  symbols  are  the  same  as  in  manuscripts 
[2,  3] . 
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W,  9jvtr )  _  WW  dm,  _  9wx  \  w  8/  ~  Bt 

r  9x1  V7n'9xi  9x1  '  9x1  9x1  ~ 

d(V*r)  _  W'  9(V,r)  W*  B(Vtr)  _  ■  _  _ 
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(+) 

O) 


For  perpendicular  intersection  (y,  z,  $)  coordinate  system  used  in  Fig.  1, 
there  are 


and 


•\/«u  *  V«a  *  1»  $u  —  90°,  Wx  —  «»,  —  IFy, 

IF*  —  —  W„  /,  —  f,  —  F„ 

/  -  a  +  [(*0*  -  («*r)*)/2,  B  -  r/r  ,  Til  -  iA  -  dp/p 


After  using  the  various  relationship  equations  as  above,  the  continuity 
equation  and  motion  equation  are,  respectively: 
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In  equations , 
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Using  p«pRT  and  W  =W  tan  6,  we  derive  from  eq.  (2*)  and  eq.  (4'): 
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Fig.  2.  Concerning  3r/3y  and  Fig.  3.  Concerning  3y/3m  and  3z/3m. 

dr  Id*  '• 


We  know  from  Fig.  2: 
We  know  from  Fig.  3: 


dr /Ox  ”  tin  E  ,  Or/dy  «■  cos  e , 
0>/dm  —sin*,  ds/d*  -  cot  5, 


Besides, 

+  dung  d« 
dy  d>»  dt  dm 


By  using  the  equations  above,  we  obtain 


Oy  dz  \  4,*)* 


(8) 


We  substitute  eq.  (7)  into  eq.  (3')  and  use  the  above  equation  (eq.  (8));  we 
obtain 


1  Op 

P  Oy 
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Equation  (9)  is  the  principal  equation  we  require. 
If  we  let 

ua4,  -  F,/Fft  ua  ft,  -  FJF9 
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and  use  eq.  (5) ,  we  obtain  Ff  -  t.oi,  .  d^V^/rdt,  F,  -  unM,  ■  d(V,r)/rdt, 
Substituting  the  above  equation  into  eq.  (9),  we  obtain; 
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For  the  computation  station  in  the  gap  of  a  blade  array,  there  are 

dB/dt  =  0  and  d(VQr)/dt  =  0 

After  substituting  the  above  equation  into  eq.  (9),  we  obtain 

Iff"  ~  -~(***+  t.n»co.e) 

+  1(1  -  M^cot£  +  Wiun^dn* ] 

(10) 

\-M\c,  dt 

Neglecting  the  last  term  of  eq.  (10),  this  is  eq.  (1).  We  know  from 
this  that  the  equation  used  by  Heterington  is  correct. 

TRANSFORMATION  OF  Y  DIRECTION  MOTION  EQUATION 

Through  transformation  of  eq.  (9),  we  can  obtain  eq.  (11)  of  [2]. 

This  explains  that  the  difference  in  form  between  the  two  is  due  to  different 
coordinate  systems.  Actually,  these  two  are  in  complete  agreement. 


We  know  from  Fig.  1 
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Fig.  4.  Concerning  3y/3r  and  3z/3r. 


From  Fig.  4, 
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We  obtain 

dtantf  _  coi*^  dun  8 
dy  co**tr  dy 


Substituting  the  above  equation  into  the  expression  equation  of  3tano/3r, 
obtain  the  following  equation  after  putting  into  order: 


dtan8  M  CMff  d  tan  O'  .  coi*g/  ^*r\ 

dy  cos  8  dr  cM*8Wa*'# 


05) 


For  the  blade  force: 


F,  —  F,  cote  —  F.aae 
F,  —  F,one  +  F,  cose 

If  we  let  tan i  -  F,/F„  tan p- FJF„ 

By  using  eq.  (5) ,  we  obtain: 


(16) 

(17) 


Ff-  tani  •  J(V,r)/(rdt)  (18) 

F,  «■  tan/i  •  d(V9r)/(rit)  (19) 

Substituting  eqs.  (11)  through  (19)  into  eq.  (9),  we  obtain  eq.  (11)  of 
[2].  Hence,  eq.  (9)  is  in  agreement  with  eq.  (11)  of  [2]. 


The  author  expresses  gratitude  to  colleagues  Li  Genzhen,  Chen  Naixing, 
Qiang  Guofang  and  Hong  Hejing  for  their  considerable  contribution  during 
discussions . 
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APPENDIX  DERIVATION  OF  TRANSFORMATION  OF  EQ.  (9)  INTO  EQ.  (11)  of  [2] 


We  obtain  from  eqs.  (11)  through  (13), 

lias  +  tin  5  cm*  —  iiaff/cm  t  (A-l) 

(1  —  IT,)  cot  1  +  VJ  naff  tie*  —  (l—  VJ)eo»»  —  VJ  naffiia*  (A-2) 

tfj  Hteeil  —  (1  —  M})iiat  ”  VJ  tmnO cm t  4-  (taa'aJWJ  —  l)ua«  (A-J) 

Substituting  eqs.  (14) - (17) ,  (A-l),  (A-2)  and  (A-3)  into  eq.  (9),  we 
obtain  after  ordering: 
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And 
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Substituting  eqs.  (A-5)  to  (A-7)  into  eq.  (A-4)  and  considering  eqs. 
(12),  (18)  and  (19),  we  obtain 
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Neglecting  the  term  ds/da  in  eq.  (A-8)  and  considering 


(  g  tan  g  \ 

V 7?7#  1 

^  *  K 

We  obtain  eq.  (11)  of  [2].  Therefore,  these  two  equations  are  in  complete 
agreement . 


(11), 
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NUMERICAL  CALCULATION  OF  SUBSONIC  COMPRESSIBLE  TURBULENT  FLOW  IN  A  CONICAL 
DIFFUSER* 


Gao  Lijun,  Wu  Wendong,  and  Wang  Yingshi 

Institute  of  Engineering  Thermophysics,  Chinese  Academy  of  Sciences 


In  this  paper,  vorticity  <a  and  stream  function  V  are  used  as  dependent  variables 
to  predict  the  turbulent  flow.  Generally,  the  method  can  be  used  to  calculate  both 
recirculating  flow  and  non-recirculating  flow.  The  calculation  of  the  turbulent  flow 
in  a  conical  diffuser  having  a  total  divergence  angle  of  10°,  and  with  on  uniform 
velocity  distribution  at  entry  is  carried  out,  Because  the  flow  relates  only  to  the 
developing  process  of  turbulence  boundary  layer,  Prandte  mixing  length  hypothesis 
is  used  in  this  calculation. 

The  numerical  results  agree  with  the  experimental  data'"  quite  well.  It  shows 
that  this  method  can  be  applied  to  the  design  of  diffuser  without  recirculating  flow. 


I .  FOREWORD 

For  more  than  ten  years,  the  numerical  calculation  method  was  applied 
in  the  flow  process  of  a  diffuser  because  of  rapid  development  of  calculation 
fluid  mechanics.  The  paper  tries  to  begin  from  an  elliptical  model  equation 
to  solve  for  the  viscous  flow  problem  in  a  diffuser.  Theoretically,  this 
approach  can  be  adapted  to  solve  for  a  flow  situation  with  or  without 


♦The  paper  was  read  at  the  Third  All  China  Engineering  Thermophysics  Confer¬ 
ence  at  Guilin  in  April  1980.  33 


recirculating  flow.  As  the  first  step,  the  authors  calculate  the  flow  of  a 
small  expansion  angle  diffuser.  Since  the  flow  is  a  development  process  of 
a  turbulent  flow  boundary  layer,  the  calculation  of  the  turbulent  flow 
viscosity  coefficient  can  use  the  Prandtl  mixing  length  hypothesis.  The 
expansion  angle  of  the  diffuser  is  10°,  and  the  data  used  for  comparison  are 
taken  from  manuscript  [4] .  The  calculation  result  and  experimental  data  are 
in  good  agreement. 

II.  FUNDAMENTAL  EQUATIONS 

The  authors  use  the  spherical  coordinate  system  (51=R,  ?2=6,  and  ^2=e^’ 
referring  to  Fig.  1.  The  fundamental  equation  set  is  established  by  focusing 
on  stability,  axial  symmetry  and  two  dimensional  flow.  In  addition,  assump¬ 
tions  include  that  the  ideal  gas,  and  C  and  A  are  constants. 

P  . 


Fig.  1 

The  continuity  equation: 

Q(RrGt)/dR  +  d(rGi)/d(l  —  0  (1) 

Gj  and  are  components  of  G=pV. 

The  momentum  equation: 

dw(CVj  —  TO  —  (fiVl  -  TvO/R  +  (Tu/fjdn^  +  Op/dR  -  0  (2) 

ft:  -  TO]  -  (Tht/r)e otp  -  (GtV,  -  TM)  •  (dR/dfi)/R> 

+  (dp/dfi)/R  -  0  (3) 

In  order  to  avoid  the  difficulty  of  determining  the  pressure  boundary 
conditions,  vorticity  w  is  introduced  in  the  calculation,  and  the  pressure 
terms  in  eqs.  (2)  and  (3)  are  cancelled,  in  deriving  the  momentum  eq.  [1] 
with  w  as  the  variable.  34 


^{a[(*./r) .  (S4/a*)]/a*  -  a[(«/r)  •  ( e*/a*)]/a* ) 

-  aUr'dt^Oa/OJ/aaj/a*  -  d{(r>/R)  •  a^o/oya^/a* 

-  r*{dL(Vl  +  VD/2]/dR}(dp/dfi) 

+  ^{a[(F{  +  Vi)/2]/dfi}(dp/dR)  -  0  (4) 

The  flow  function  equation: 

a{[*/(pOKa*/d/?)}/a/?  +  d{[i/(P*,)Kd*/a/0}/a/?  +  ru  -  o  (5) 
The  energy  equation: 


e(*»a*/8£)/aR  -  a(A*a</>/a*)A»  -  d(jk,MRrQh*/dR)/dR 

-  d[Tk*tt(r/R')dh*/ae)/dp  -  dfp^KrO  -  Uak^)d{Vii2)/dR\/dR 

-  a{[M^i(r//J)  •  (1  -  \/<J>«*)]Q(V'/2)/de)/dfi  -  0  (6) 


In  eq.  (6),  h*  is  the  total  enthalpy;  e£.  is  the  effective  Prandtl  number 
of  the  turbulent  flow;  is  the  effective  heat  exchange  coefficient  of 

the  turbulent  flow;  and  ^  eff=Meff/ah  eff>  ye££  is  the  effective  viscosity 
coefficient  of  the  turbulent  flow. 


From  eqs.  (4)  to  (6),  the  general  form  of  the  fundamental  equation  can 
be  derived: 

«♦£ OC+dt/dri/dR  -  d(4>d<l>/QR)/dfi]  -  3[bJird(c.+)/dR}/dR 

-  a[A*(r/*)d(^)/a/»]/a^  +  Rrd.  -  o  (7) 

Corresponding  to  different  coefficients  a^,  b^,  c^,  and  d^  of  different 
variables  <t>,  the  following  table  shows  the  contents: 


4> 

•* 

*♦ 

<• 

to/r 

r' 

r* 

Ptii 

-  ('•/R){<«p/a/»>5Kn  +  K])/2J/<Mt  -  (ep/a*)«t(K}  +  v\)/2]iap} 

4> 

0 

KpH) 

1 

—  (u/r 

*• 

1 

1 

-  n/R)d[»M(Mi  -  1 /«..*.)«< v'/2)/<«] /a* 

-  0/^)atp.fl(r/RXi  -  i/»*,wi)a(t"/2 )ldPl/»fi 

For  solving  variables  w/r,  i>  and  h*,  the  calculation  of  the  related  physical 
quantities  is  as  follows: 
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Velocity: 


V\  -  (d*/9fi)/(pRr),  V>  -  -(8*/d*)/ (pr) 


(8) 


Pressure:  Calculation  of  pressure  is  based  on  the  momentum  equations 

(2)  and  (3);  the  calculation  of  the  relative  viscosity  stress  in  eqs.  (2)  and 

(3)  is  as  follows: 


Tut  -  f^aUdVjdR  -  (2/3)divVH  (9) 

f^[2 {dVjdfi)/R  +  IV JR  —  (2/3 )div  V]  (10)’ 

Tm  “  ^«M[2(Fisin^  ■+■  V2co»0)/r  —  (2/3)divV]  (11) 

Tu i  -  Tj.,  -  f^[{dVj6fi)/R  +  Rd{VjR)/dR]  (12) 


In  the  equations, 

divV-  1  /{Rr)[RrdVx/dR  +  FiRsifl*  +  Vxr  +  rdVjdfi  +  V}Rco*0]  (13) 

We  substitute  eqs  (9)  to  (13)  into  the  momentum  equations  (2)  and  (3); 
3p/ 3R  and  3p/36  can  be  derived  with  values  shown  in  the  following: 


dp/dR  ■=  -  (l/i?r)d{/?r[^  -  l^QVjdR  -t-  (2/3)/Urff  div  V]  }/dR 

-  (1  /RrXdlrpVS,  -  r^{dVJQfi)/R  -  r^QVjQR 
~  VJR))/Qfi)  +  [PV\  -  l^dVJdft/R 

-  Ip^VjR  +  (2/3  WivV]//? 

-  [2f*»«(  Firing  -+-  F2cos  P)/r 

~  (2/3)^,  div  V)/R  (14) 

d/»/a/?=  [i/(i?r)]d{/?2r[^(aF2/ai?  -  vjr)  +  m(dvjdp)/R) 

—  pR2rVlVi}/dR  +  U/{Rr)}d{Rr[^t{{2/R)dV1/dfi 
+  lp*tV\/R  —  (2/3)^  div  V  —  pV\]  }/Qp 

—  f*.«[2( Firing  +  F2co *fi)/r 

—  (2/3)divV]/un0  (15) 

P*  -  Pa  +  0.5 [{dp/dR)A  +  (dp/dR ),](**  -  *,)  (16) 


or  pc  —  po  +  0.5[(a?/d/»)e  •+■  (dp/dfi)DKfic  ~  to) 

Density:  P~p/gRT 

In  the  equation,  T  •*  (hm  —  V2/2g)/c, 

Effective  viscosity  coefficient  of  turbulent  flow: 
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07) 

(18) 

(19) 

Pi  + 


(20) 


f  l 


is  the  viscosity  coefficient  of  the  laminar  flow;  and  Up  is  the  viscosity 
coefficient  of  the  turbulent  flow.  The  calculation  of  u  adopts  the  Prandtl 
mixing  length  hypothesis ;  i.e.,  Mr  “  »/®>  !• 

The  calculation  of  mixing  length  1  is  based  on  Escudier's  empirical 
When  y/5$  A/k,  m 

,  /. .  ( 21 )  formula  5  . 

When  y/6>  A/k,  !./•-*>/>  ( 

ija  - 1  (22) 

In  the  equations,  k=0.41  and  X=0.09. 

In  the  two  above  equations,  y  is  the  vertical  distance  from  the  cone 

surface  of  the  diffuser;  6  is  the  boundary  layer  thickness.  In  the  calculation, 

V./V  >0.99  is  used  to  determine  6.  From  eqs.  (21)  and  (22),  it  is 
i  max 

apparent  that  there  is  a  certain  limitation  in  using  the  Prandtl  mixing  length 
hypothesis.  When  \&V\/Qy |  -»■  0, 

lip  seems  to  have  the  tendency  of  approaching  zero;  of  course,  this  is  not 

consistent  with  the  physical  significance,  so  this  is  treated  in  calculation 

[2].  When  |3Vp/3y|  is  smaller  than  a  certain  small  quantity  (corresponding 

to  V./V  5.0.99)  and  if  we  let  /«*  **  plm(Ll8Vi/dy\)- 

i  max 

In  the  equation,  /.|a*yayj  -  0.01(FW  +  Vui«)/2, 


j  and  j+1  indicate  two  adjacent  nodal  points  along  the  direction. 


III.  BOUNDARY  CONDITIONS 


Inlet:  co/r  -  (RfiVJdR  -  The  lower  angle  0  indicates  the 

inlet . 


V  “  Rlfit  ( 1  P1(oG?)sin£</|J,  V,  n  is  the  non-uniform  distribution; 

taken  from  [4],  V.,  g=0. 

**  “  CfT0+  V\*l2g,  293K  is  taken  for  T  . 

Upper  wall  surface:  (* /r )y  -  _  [3(V.v_,  -  0  +  0.5(a»/r)v_t) - 

The  lower  angles  N  and  N-l  indicate  the  nodal  points  on  the  wall  surface,  and 
the  nodal  point  of  a  mesh  lattice  separating  from  the  wall  surface. 


dhm/dy  —  (VfiVj/dy^/g  —  0,  y  is  the  vertical  distance  from  the  wall  surface. 
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The  central  line:  „/,  -  -  8[(*,  -  -  (9t  -  *,)  ^]/[PR\fii  -  #)}. 


The  lower  angles  1,  2  and  3  indicate  the  central  line  and  the  adjacent  nodal 
points,  referring  to  Fig.  2.  -  0; 

0utlet  :  d(ta/r)/dR  —  QV/dR  —  dh*/QR  —  0 


h2 

-  J - ***(*) 


Vc 


Fig-  2.  Fig.  3. 

Key:  (*)  Central  line.  Key:  (a)  Line  of  equal  8;  (b)  Curve  of 

equal  R. 


IV.  ESTABLISHING  AND  SOLVING  FOR  DIFFERENCE  EQUATIONS 


Derived  from  eq.  (7),  the  difference  equation  solving  for  numerical 
calculation  should  be  able  to  plot  a  series  of  mutually  perpendicular  lines 
of  equal  R  and  lines  of  equal  6  along  directions  of  R  and  8.  We  integrate 
the  control  area  enveloping  the  point  F  (the  portion  of  cross  sectional 
lines  in  Fig.  3).  In  addition,  in  the  convection  term  integration,  the 
average  variables  are  taken  as  the  up-wind  value  in  deriving  the  following 
difference  equation  [1]: 


y,  {[A]  + 

LzHiLLS. - 

y  {A]  ■+■  ?♦,/•(**„  ■+■ 

/•.v.j.i.ir 


(23) 


In  the  equation, 

A ;  -  A,/V„  B ;  -  B,/lVf(Ki  +  *„,)]: 

At  -  («*,/»)l (*«  +  V,  -  Vst  ~  +  |*«  +  Vs  -  r.vs  -  ^.vl  ] 

-  (««rf/S)[(Vwr  +  Vs  -  Vs w  -  vs). A-  | Vsw  +  Vs  -Vsw-Vs\] 
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as  -  (#*.,/« )[(y,»  +  r.  -  -  *,)  +  |vfr  +  »»  -  »,£  -  r,|  ] 

•^N  “  (<**.f/8  )l(V\-r  +  &£  —  9s»-  —  $rM  )  4-  ]r.v,  +  Vt  —  Vyr  —  •Vj] 
Bf  “  Rr(b*. e  4-  ^t.r'Kfis  £j)(Br  +  B?)  sin  ^f/[8(Bf  B^)J 
B»  "■  Rr(b+,a-  4-  b+.r'Xfs  ~  Ps)(Rw  4-  B*)sin£*/(8(B*  —  B*/)} 

By  “  (i*.v  4-  i*,f)(Br  —  B»X nn Bn  4-  sin  £?)/[8(£v  —  |ff)] 

Bs  *=  (*♦.*  4-  A*,fXBr  Bg  )(  sin  4-  sin  Pr)/[%CPp  £s)l 

F,  -  R\aap,(Rt  -  RwXfix  ~  ?$)/< 


For  convenience  of  convergence,  the  implicit  iterated  substitution 
formula  is  taken  for  the  w/r  value  on  N-l  points. 

y,  (  [Aj  4"  4-  i«,.v-i)By}(^)( }  4-  {  }yCi —  ^*,y-i 

(«/0v-,  -  i=^— - - - - - ,  (24) 

/  .  4-  <"*,y-i( £«,/  4-  ^,v-t)B/}  —  {  }yC2 

l*s.s.t.w 

<}>  represents  the  corresponding  cj/t. 

{  }v  *  {<<*■  4  <♦,!*($*,*  4-  A*.*)B'„.} 

~  ~  Xf.v-:  ”  )/ ( ^y-iny-tPy-! ) »  «*_,  —  B(£y  — 

C2  -  —  1/2 

Selections  of  mesh  lattices  are  as  follows: 

Sj-.  equidistant  mesh  lattice  with  lattice  distance  of  0.02  meter. 

^2:  unequal  mesh  lattice.  For  viscous  flow,  since  a  development  process 
exists  in  the  boundary  layer,  the  velocity  variation  is  greater  within  the 
boundary  layer,  so  the  mesh  lattices  should  be  denser.  The  mesh  lattices 
which  are  wider  apart  will  cause  inaccurate  calculation  results. 

Figure  5  indicates  the  influence  on  the  central  velocity  of 
[e  —  (£.*,  —  —  fii-0  —  0.9,  0.943,  1  ]  for  selecting  different  values  of  e  . 

Compared  to  the  experimental  values,  the  error  is  the  smallest  whene*0.9. 
Therefore,  the  calculation  proceeds  with  e*0.9. 
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In  this  calculation,  the  superrelaxation  method  is  used  for  solution; 
the  superrelaxation  factors  are:  a„,«*  1,  a j  *  1.25,  a*»  —  l  Calculating 

for  convergence  on  all  internal  points  satisfying 

The  value  of  X  is  0.0005;  n  represents  the  number  of  repetitions  at  this 
time;  and  (n-1)  represents  the  number  of  repetitions  for  the  last  time. 

V.  CALCULATION  RESULTS  AND  ANALYSIS 

This  calculation  is  conducted  on  a  straight  wall  conical  diffuser  with 
divergence  angle  of  10°,  inlet  radius  of  0.07762  meter,  and  length  of  0.2  meter 
proceeding  under  the  condition  of  uneven  velocity  at  the  inlet.  The  calcula¬ 
tion  results  are  as  follows: 

Figure  4  indicates  the  variation  of  the  static  pressure  restoration 
coefficient  C?  «  (pi  —  fl»)/(0.5pflV?,fc.)  with  R.  The  static  pressure 

is  the  arithmetic  mean  of  pressures  at  various  points  of  the  cross  sectional 
surface.  q  m  *n<*icates  the  fringe  velocity  at  the  boundary  layer  of  the 
inlet  along  the  ^  direction. 

Figure  5  indicates  the  situation  that  the  velocity  at  the  fringe  of 
the  boundary  layer  reduces  with  increase  in  diffuser  length.  Compared  to 
the  experimental  data  in  [4],  this  is  relatively  satisfactory. 

Figure  6  indicates  the  influence  on  value  5/Rq  of  two  inlet  velocities 
(uniform  and  non-uniform).  The  calculation  is  conducted  in  the  situation 
that  the  inlet  volumetric  flow  is  under  the  same  conditions.  The  inlet  non- 
uniform  velocity  distribution  (Fig.  7)  is  derived  using  the  method  provided 
by  [3]  with  the  fringe  velocity  of  the  boundary  layer  as  provided  by  [4] . 

From  Fig.  6,  it  is  apparent  that  the  difference  is  relatively  high  when 
induced  by  two  different  inlet  velocities. 


Figure  9  indicates  the  situation  when  the  turbulent  flow  viscosity 
coefficient  at  a  certain  cross  sectional  surface  varies  with  6.  In  the 
diagram,  values  of  points  A,  B  and  C  are  apparently  irrational;  this  is 


caused  by  using  the  Prandtl  mixing 
tion  result,  there  is  no  undesirabl 
(Fig.  8). 


length  model.  However,  from  the  calcula- 
e  influence  on  the  velocity  distribution 


l'*‘  0.»a  0.92  0.94  1.00  1.04  1.08 

*(*)  (a) 

A  o  ▼  if  * 

«: 

Fig.  4.  (e) 

Key:  (a)  Meter;  (b)  Experimental 
value;  (c)  Calculation  value;  (d) 
Calculation  value;  (e)  p  takes 
part  in  repetitive  substitutions 
intermittently. 


u  97  >  *ft«(b)  04o  O  o  4 

0  5 ! _ i--  l»CJ  ■  . 

0.88  0.92  0J4  18  l.#4  i.o*  i.j# 

-  *C*)(a) 

Fig.  S. 

Key:  (a)  Meter;  (b)  Experimental  value; 

(c)  Calculation  value. 


*  *  ***  *£»«>**  «/» 


041  042  0.94  1.00  1.04  1.01  1.10 

a  o  tfici  * 

x  tf*. 

ft:  *S*a**a>*,**#w 

Fig.  6.  .  .  *#*«,# 

Key:  (a)  Meter;  (b)  Experimental 
value;  (c)  Calculation  value: 
velocity  distribution  of  nonuniform 
inlet  velocity;  (d)  Calculation 
value:  distribution  of  uniform 
inlet  velocity,  distributed  accord¬ 
ing  to  equal  volume  flow. 


0480  0.085  0.0*0 

KXm'i  (a) 


Fig.  7. 

Key:  (a)  Radian;  (b)  Velocity  of  equal 
volume  flow  uniformly  distributed. 


VI.  CONCLUSIONS 


(1)  If  we  compare  the  calculation  results  using  this  method  and  that  of 
the  experimental  results,  the  degree^f  agreement  is  good. 


Fig.  8.  Fig.  9. 

Key:  (a)  Calculated  value;  (b) 

Manuscript . 


(2)  The  calculation  proves  the  important  influence  on  flow  field  of  the 
inlet  velocity  distribution.  The  calculation  results  are  obviously  different 
using  uniform  and  nonuniform  inlet  velocity  distributions  on  the  same  diffuser. 


(3)  This  program  can  be  used  to  calculate  the  flow  in  a  small  expansion 
angle  diffuser  of  the  same  category  flow  without  separation.  As  the  compressi 
bility  influence  is  counted,  this  program  can  be  used  for  calculation  of 
subsonic  flow  of  a  wider  range  of  inlet  M  number. 


(4)  To  shorten  the  calculation  time,  in  addition  to  the  adoption  of 
better  relaxation  factors,  the  pressure  subprogram  can  intermittently  take 
part  in  repetitive  substitutions.  Thus,  the  accuracy  of  calculation  results 
is  not  affected,  referring  to  Fig.  4. 


The  authors  express  gratitude  to  Professor  Wu  Zhonghua's  guidance. 
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THE  ANALYSIS  OF  FORCES  AND  POWERS  IN  TURBOMACHINERY  GAS  DYNAMICS* 
Wang  Jiasheng 

Shanghai  Institute  of  Mechanical  Engineering 


Correctly  analyse  the  forces  and  Powers  are  an  very  important  problem  in 
tnrbomachinery  gas  dynamics.  This  problem  was  discussed  detaiiy  in  this  paper  and 
was  concluded  that  the  realization  of  forces  and  powers  are  concerned  with  the 
situation  of  observers.  When  we  disease  the  forces,  we  must  point  ont  the  situation 
of  observers ;  When  we  discuss  the  powers,  we  must  point  out  what  force  to  work  in  . 
what*  motion. 

The  paper  is  written  on  the  basis  of  the  characteristics  of  gas  motion 
in  turbomachinery  and  the  situation  of  observers;  the  paper  stresses  the 
discussion  of  acting  forces  and  the  realization  of  power  produced. 

We  arbitrarily  take  a  micro-element  control  body  dV=dyjdy7dy3  in  a 
Cartesian  coordinate  system  oy^y^yj.  ^  i-s  no^  difficult  to  apply  the  momentum 
principle  to  derive  the  micro-element  power  dL  of  an  acting  force  on  micro¬ 
element  system  pdV,  which  is  dL  _  JL  (C.^iV  +  pfiCdV 

dy, 


We  develop  the  above  equation  and  apply  a  N-S  equation;  we  can  immediately 

derive  the  micro-element  power  of  an  acting  force  on  a  unit  of  mass  of  gas 

dl  -  dL/pdV  -  C-  ~  +  fnDt,  O) 

Dl  P 


•The  paper  was  read  at  the  Third  All  China  Engineering  Thermophysics  Conference 
at  Guilin  in  April  1980.  43 


In  the  equation,  D^=3Cj/3y^  is  the  element  of  the  displacement  tensor. 


It  is  apparent  from  eq.  (1)  that  the  micro  power  (acting  force  on  a  unit 
mass  of  gas)  is  composed  of  two  parts:  the  first  is  the  numerical  product 
of  velocity  C  and  acceleration  DC/Dt;  this  is  not  related  to  deformation  of 
the  gas  micro-element  system.  This  portion  of  the  power  is  called  the 
migrating  power.  The  second  is  related  to  the  stress  tensor  and  displacement 
tensor;  i.e.,  it  is  related  to  deformation  of  the  gas  micro-element  system. 
This  portion  of  the  power  is  called  the  deformation  power. 


For  the  migrating  power,  since  the  unit  mass  of  gas  is  taken,  the  force 
acting  on  a  unit  mass  of  gas  is  equivalent  to  acceleration  of  gas  motion. 
Then,  the  discussion  of  the  acting  force  can  be  substituted  with  discussion 
of  acceleration. 


As  is  well  known,  the  realization  of  velocity  and  acceleration  in  turbo¬ 
machinery  is  related  to  the  observer's  position.  Or,  it  is  necessary  to 
point  out  whether  the  adopted  system  is  inertia  system  or  not  an  inertia 
system.  For  convenience  in  discussion,  we  call  the  observer  situated  in  an 
inertia  system  the  absolute  observer,  while  the  observer  situated  in  the  non¬ 
inertia  system  is  called  the  relative  observer. 

In  the  following,  the  realization  problem  of  the  acting  force  is  discussed 
first. 

(A)  The  problem  of  the  acting  force  of  a  gas 

(1)  To  an  absolute  observer,  the  observed  absolute  acceleration  D  C/Dt 
can  be  expressed  as 

DjC/Dt  -  DW/Dt  +  2*  X  W+  0  X  (»  X  r,)  (2) 

It  is  apparent  from  the  above  equation,  as  viewed  by  an  absolute  observer, 
that  the  force  acting  on  a  unit  mass  of  gas  include:  (a)  the  force  DW/Dt 
producing  relative  acceleration  of  a  unit  mass  of  gas,  (b)  the  force  2«»x  W, 
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producing  Geshi  [transliteration]  of  a  unit  mass  of  gas  is  called  the  Geshi 
force;  and  (c)  the  force  wx(a>xrp)  producing  centripetal  acceleration  of  a 
unit  mass  of  gas  is  called  the  centripetal  force. 

(2)  To  a  relative  observer  moving  with  the  rotating  wheel,  the  observer 
views  the  rotating  wheel  as  stationary  but  the  gas  (in  the  rotating  wheel) 
moves  with  velocity  W  and  acceleration  DW/Dt.  However,  the  force  the  observer 
actually  observes  is  D  C/Dt,  D  C/Dt^DW/Dt, 

ci  cL 

that  is,  the  force  acting  on  a  unit  mass  of  gas  is  not  equivalent  to  the 
acceleration  thus  produced.  In  order  to  understand  this  conflicting  phenomenon, 
the  observer  assumes  two  other  forces  as  he  considers  that  the  acting  force 
should  include:  (a)  the  force  D  C/Dt  directly  observed  by  the  relative  and 

cL 

absolute  observers;  (b)  the  assumed  inertia  force  -wx(wxrp),  the  so-called 
centrifugal  force;  and  (c)  another  assumed  inertia  force  -2wxW.  In  order  to 
distinguish  the  force  from  the  Geshi  force  in  the  inertia  system,  this  force 
is  called  the  Geshi  inertia  force.  Then 

DW/Dt=DaC/Dt/w  x  (w  XTp)  -2w  xW 

It  is  apparent  that  while  analyzing  the  action  force  in  the  gas  dynamics 

of  turbomachinery,  it  is  necessary  to  mention  the  position  where  the  observer 

is  attempted.  To  an  absolute  observer,  the  Geshi  force  and  centripetal  force 

make  sense  because  these  forces  are  observed  by  the  absolute  observer  as  two 

force  components  of  a  force  D  C/Dt.  To  a  relative  observer,  the  Geshi  inertia 

force  and  centrifugal  force  make  sense  because  these  two  forces  explain  to 

the  relative  observer  an  assumed  force  under  the  action  of  a  force  D  C/Dt  but 

a 

have  an  "irrational"  phenomenon  of  producing  DW/Dt  acceleration. 

Together  with  this  conclusion,  there  necessarily  exists  a  problem  of 
views  on  producing  power  by  two  observers. 

(B)  The  problem  of  producing  power  by  gas 

(1)  Migrating  power  From  equation  (1)  and  C*W+u>xrp,  we  can  obtain 

C-D  C/Dt« [DW/Dt  +<iJx(wxr  )  +  2uxW]  •  (W-^xr  ). 
a  P  P 
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The  left  side  of  the  above  equation  can  be  rewritten  into  (D  /Dt)  (C2/2) .  On 
the  right  side,  six  terms  can  be  obtained  after  development.  Their  individual 
physical  significance  can  be  discussed  as  follows: 

indicates  the  power  created  by  the 

force  producing  the  relative  acceleration  for  the  relative  motion.  The  term 
is  equivalent  to  the  variation  rate  with  time  t  of  relative  kinetic  energy 
by  a  unit  mass  of  gas. 

(b)  DW/Dt  •  (<p  X  r,)  —  m  indicates  the  power 

produced  by  force  DW/Dt  for  tractive  motion. 

(c)  [»  x  (a»  x  r,)]  •  W  —  — —  indicates  the  power  acting  on 

the  gas  by  a  centripetal  force  in  relative  motion. 

(d)  2  (oixW)  •  W=0  indicates  the  power  acting  on  the  gas  by  the  Geshi  force 
in  relative  motion.  Since  the  Geshi  force  is  consistently  perpendicular  to 
the  relative  velocity  W,  the  power  produced  is  zero. 

(e)  2(«  x  HQ  •  (•  X  r()  ■  2  —  (-)  indicates  the  power  acting 

Dt  \  2  / 

on  the  gas  by  the  Geshi  force  in  tractive  motion;  this  indicates  that  the 
Geshi  force  does  produce  power  on  tractive  motion. 

(f)  [u)x(u«rp]  ]  •  (wxrp)=0  is  the  power  produced  (zero)  by  centripetal  force 
on  tractive  power. 

Similar  to  discussion  of  the  acting  force,  to  a  relative  observer, 
since  C*D  C/Dt^DW/Dt-W,  in  order  to  explain  this  ’’irrational”  phenomenon,  in 
the  view  of  the  relative  observer,  in  addition  to  DW/Dt ‘W.  the  virtual  power 
should  be  added:  -DW/Dt • (wxrp) ;  [-wx(u)xrp)]  ’W;  [-2u)xW]*W. 


Comprehending  the  above  analyses,  we  come  to  the  view: 


(a)  While  discussing  the  power  created  on  the  gas,  it  is  necessary  to 
make  known  the  position  the  observer  is  in.  To  an  absolute  observer,  the 
Geshi  force  and  centripetal  force  make  sense  for  power  acting  on  the  gas. 

To  a  relative  observer,  the  Geshi  inertia  force  and  centrifugal  force  make 
sense. 

(b)  It  is  irrational  to  make  the  statement,  "a  certain  force  does  work 
(not  does  'no'  work)"  but  to  clarify,  "a  certain  force  does  work  (or  does  no 
work)  on  certain  motion." 


(2)  The  deformed  power  can  use  the  stress  tensor  expression  equation 
of  equation  (1)  and  Newtonian  fluid  pu  —  —  p9it  +  2^  —  2/3pS„  div  C 


rewritten  into  the  following  (in  the  equation 
S_  is  the  deformation  tensor  element) : 


(S  .  . 


is  the  Kronecker  sign  and 


From 

We  obtain 


*■  O/pX”  p9»  +  2 fiSa  —  2/lftdtj  div  C)D# 
Su  ”  Sin  A#  ”  On,  Da  ”  Sj,  +  Qa 

9nD k  —  Du  ™  div  C,  SifOa  “  0 

(1  /APhDh  -  0/p)[— pdiv  C  +  2^5,  -  (2/3)p(div  C?] 


Let  jP  -„[2S!,-(2/3XdivO’],  then 

(l/p)ft/Av  —  —  (/>/p)dW  C  +  P/fi  (3) 

In  the  equation,  <P  is  the  divergence  function. 


For  power  produced  in  tl.is  portion,  we  first  notice  that  p,  p  ,  and 
div  C  are  not  related  to  the  position  of  the  observer.  Hence,  it  can  be 
considered  that  the  deformation  power  is  not  related  to  the  observer. 
Secondly,  we  know  from  the  divergence  definition  formula  div  C  =*  6  (dV)/dVdt 
of  velocity  C,  -(p/p)div  C  =  -p[6  (dV)]/pdVds.  In  the  equation,  5  (dV) 
indicates  the  variation  on  micro-element  volume  dV  in  dt.  Hence,  this  term 
can  be  viewed  as  the  micro-element  compression  power  on  a  unit  mass  of 
gas  by  pressure. 
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FLOW-FIELD  LINE  RELAXATION  SOLUTION  FOR  S?  RELATIVE  STREAM  SURFACE  WITH  THE 
SPLITTER  BOUNDARY  CONDITION* 

Zhu  Rongguo 

Institute  of  Thermophysics  Engineering,  Chinese  Academy  of  Sciences 

Based  on  Wu’s  general  three-demensional  turbomachine  flow  theory,  &  computer 
program  of  flow  field  matrix  line  relaxation  waa  established  for  flow  along  the 
relative  stream  surface  Si  of  ducted  fans  in  bypaaa  engine. 

It  was  the  development  of  reference  (3),  and  the  splitter  boundary  condition  of 
the  fan  compressor  was  considered  and  the  structure  of  the  computer  program  was 
improved. 

I.  FOREWORD 

In  a  fanjet  engine,  the  shape  and  position  of  splitters  of  internal  and 
external  inclusion  can  have  a  significant  influence  on  performance  of  a 
fanjet  engine.  Hence,  the  understanding  of  the  distribution  of  velocity 
and  pressure  of  the  splitter  surface  and  the  influence  on  the  flow  field 
have  very  important  significance  to  raise  the  performance  of  a  fanjet  engine. 

The  computer  program  introduced  in  the  paper  adopts  the  fundamental 
equation  set  [1,  2]  of  turbomachinery  three-dimensional  flow  gas  dynamics 
with  arbitrarily  non-perpendicularly  intersecting  curve  coordinates  and 
the  corresponding  non-perpendicularly  intersecting  velocity  component,  as 
well  as  the  flow  field  line  relaxation  numerical  calculation  method  [3]. 

This  is  a  development  of  the  original  S2  stream  surface  counter  problem 

♦The  paper  was  read  at  the  Third  All  China  Engineering  Thermophysics  Confer 
ence  at  Guilin  in  April  1980.  49 


computer  program  [3].  The  program  considers  the  boundary  condition  of  a 
splitter  with  internal  and  external  inclusions;  it  is  primarily  used  in  gas 
dynamics  design  computation  of  the  fan  portion  of  the  fanjet  engine. 

Compared  with  the  original  stream  surface  counter  problem  computer 
program,  this  program  makes  some  effort  in  saving  the  internal  memory 
storage  capacity,  thus  computation  nodal  points  are  used,  using  only  the 
internal  memory  from  the  originally  allowed  some  700  nodal  points  to  some  1000 
points  without  reducing  the  computation  speed. 

II.  COMPUTATION  EQUATIONS 

(A)  The  principal  equation  to  be  solved:  The  principal  equation  to  be 
solved  is  similar  to  that  described  in  manuscript  [3];  however,  for  saving 
internal  memory  capacity  during  computation,  the  equation  is  written  into 
the  following  form: 

w  [(v>  +  via (0,  -  e2))  VZl  -  ±  f( *" sb (e,  -  6,)  +  vo  yfc] 

Var  Q(Vtr)  a H  dt  1 

“  v  l  r1  ~9?~+  a?  -  Td?-U\ 

We  introduce  the  flow  function  variable  ,  definition  dV/Qx1  —  fy/T^pV' 

and  dV/dx'  —  —  fjTnfiV1  .  In  the  equations, 

f  -  2 «r(P  -  OcoiCG,  -  0;)/(C  •  P  •  Sc). 

Then  the  principal  equation  for  solving  the  flow  function  is: 

Vfu. 
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(B)  Calculation  of  other  parameters:  The  fundamental  calculation  formula 
is  the  same  as  the  original  S?  program  [3].  Several  differences  are  explained 
in  the  following: 


(a)  Since  the  program  is  mainly  used  for  gas  dynamics  calculation  of  the 
fan  portion  of  a  fanjet  engine,  the  temperature  increase  is  relatively  small. 
Therefore,  in  the  program  it  is  specified  that  the  gas  specific  heat  k  is  a 
constant:  the  corresponding  calculation  equation  of  dimensionless  temperature 
and  density  is  T=  (H-V2/2)/H1  and  p  —  T  ~*^r  • 

In  the  equation,  and  are  constants  related  to  inlet  conditions. 

(b)  By  using  the  catching  method  solving  for  the  value  of  station  j , 
since  the  value  'I'  of  the  root  portion  of  the  external  inclusion  part  does  not 
equal  zero,  the  repeated  derivation  formula  during  calculation  is  changed 
slightly. 


In  the  process  of  cancelling  elements,  in  order  to  cancel  elements  in 
the  equation 

E 2i9rt-b  E “  B t 

and  to  rewrite  into  the  standard  form  of^^+b^^+1=gk.  For  k=2,  we  have 


then 


»,+  El,9JEl,  -  BJEl,  -  ElJJEl,. 


Correspondingly,  we  have  and  g2sB2/E22-E2^ ^/E2^. 

HI.  CALCULATION  REGION  AND  SEQUENCE  OF  INTERATIVE  SUBSTITUTION 
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The  entire  calculation  is  divided  into  two  major  regions:  fan  portion, 
internal  inclusion  portion,  and  external  inclusion  portion,  as  shown  in  Fig.  1. 

The  fan  region  is  ABCFC'B'A'A;  the  internal  inclusion  region  is 
EFG'D'C'B'E;  and  the  external  inclusion  region  is  BCDGFEB.  The  last  two  j 
stations  B'B  and  C'C  in  the  fan  portion,  respectively,  coincide  with  the 
first  two  j  stations  (B'E  and  C'F  and  EB  and  FC)  of  the  internal  and  external 
inclusions.  In  other  words,  the  BCFC’B'EB  are,  respectively,  of  two  different 
calculation  regions;  this  can  be  viewed  as  the  corresponding  points  of  two 
different  mesh  lattices.  On  these  points,  for  different  regions  the  geometric 
parameters,  such  as  V^n-.  a/«b»  —  d,')  are  different;  however, 

the  physical  quantities,  such  as  V0r,  H,  s  and  o,  are  the  same.  This  way  can 
be  viewed  as  a  simple  method  of  local  dual  mesh  lattices. 


) 


Fig.  1.  Regional  classification 
diagram  of  computation. 

Key:  (1)  Fan  portion;  (2)  Exter¬ 
nal  inclusion  portion;  (3)  Split¬ 
ter;  (4)  Internal  inclusion 
portion. 

During  calculation,  calculation  is  conducted  on  j  stations  in  sequence 
from  upper  to  lower  streams  for  each  region.  In  the  fan  region,  the  calcula¬ 
tion  is  conducted  on  the  B'EB  station  to  transmit  the  physical  quantities  of 
the  corresponding  nodal  points  of  mesh  lattices  of  line  B'B  to  the  internal 
and  external  inclusion  region.  Then,  calculation  is  conducted  in  sequence 
of  j  stations  from  C'F  and  FC  lines  of  the  internal  and  external  inclusions 
until  D'G'  and  GD  are  reached.  The  physical  quantities  calculated  from  C'F 


Fig.  2.  Meridian  streamline  distribution 
of  two  stages  of  fans  with  splitter. 

Key:  (1)  Guidance  blades;  (2)  Moving 

blades;  (3)  Static  blades;  (4)  Splitter. 


and  FC  stations  are  transmitted  to  the  corresponding  nodal  points  of  mesh 
lattices  on  line  C'C  in  the  fan  region.  Thus,  a  complete  repetitive  substi¬ 
tution  is  done  until  convergence  if  attained. 

IV.  COMPUTER  PROGRAM  AND  TRIAL  COMPUTATION 

The  computer  program  is  compiled  with  ALGOL  language  in  a  TQ-16  computer 
of  32K  internal  memory.  The  operating  process  of  the  computer  program  is 
as  follows: 

(a)  Input  the  original  data,  dimensionless  parameters,  and  geometric 
parameters  at  the  nodal  points  of  the  calculation  mesh  lattice. 

(b)  The  initial  field  for  repetitive  substitution  calculation  is  given 
before  the  repetitive  calculation  is  conducted,  including  the  transmitting 
values  between  the  fan  region  and  internal  and  external  inclusion  regions. 

(c)  Accuracy  inspection  is  done.  If  accuracy  is  not  accepted,  we 
repeat  (b)  [calculation  of  other  parameters,  in  II].  If  the  accuracy  is 
satisfied,  the  calculation  result  is  given  as  an  output  printout. 

The  source  program  using  ALGOL  language  is  not  presented  here. 

The  program  is  used  for  trial  calculation  of  some  examples.  One  of 
the  examples  is  for  the  two  stages  fan  of  JT-3D.  The  result  of  trial 
computation  satisfies  the  program. 

Figure  2  shows  the  meridian  streamline  distribution  under  the  operating 
condition  of  design  points  for  the  fan  portion  of  JT-3D.  Trial  calculations 
for  changing  the  channel  inclusion  ratio  were  conducted.  When  the  geometrical 
shape  remains  the  same  while  the  internal  inclusion  flow  is  reduced,  it  is 
obvious  that  the  streamlines  move  toward  the  blade  tip. 

During  repetitive  substitution,  the  super  relaxation  factor  a*l-l .65  is 
adopted.  For  satisfying  the  engineering  accuracy,  it  is  enough  to  conduct 


15-20  repetitive  substitutions.  The  computer  time  required  for  computation 
is  approximately  3-4  minutes. 
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EFFECT  OF  CASING  TREATMENT  ON  PERFORMANCE  OF  AN  AXIAL  COMPRESSOR* 

Li  Kerning,  Zhao  Quanchun,  and  Sun  Yuanying 
Shenyang  Aeroengine  Research  Institute 

An  experimental  investigation  waa  conducted,  and  is  herein  reported,  to  obtain 
the  effectiveness  of  four  treated  casing  configurations  on  performance  of  an  axial 
compressor.  Test  shows  that  the  surge  margin  and  distortion  tolerance  for  the  axial 
compressor  are  improved  by  means  of  casing  treatment,  so  that  the  circumferential 
pressure  profile  is  further  improved  and  the  peak  efleciency  of  the  compressor  is 
almost  unaffected.  There  is  a  common  trend  for  the  four  treated  casings,  that  is  the 
increment  in  surge  margin  at  lowers  peed  is  greater  than  that  at  higher  speed;  the  incre¬ 
ment  in  surge  margin  with  distorted  inlet  flow  is  more  than  that  with  uniform  inlet  flew. 

I .  FOREWORD 

With  increased  range  of  aircraft  flight  and  the  appearance  of  anomaly 
flow  fields  at  the  engine  inlet,  it  is  required  to  have  sufficient  distortion 
tolerance  and  surge  margin  capability  for  the  compressor.  The  casing  treatment 
is  one  of  the  effective  measures  of  increasing  distortion  tolerance  and 
surge  margin. 

*¥he"  paper  was  read  at  the  Third  All  China  Engineering  Thermophysics  Conference 
at  Guilin  in  April  1980. 

SS 


The  report  introduces  the  principal  results  of  a  three  stage  super- 
and  transonic  low  pressure  compressor  with  experimental  study  [1]  on  casing 
treatment  for  the  first  stage  rotor  outer  wall.  Four  types  of  casing 
treatment  were  designed  and  tested  (small  blade  turbulent  flow  type,  circum¬ 
ferential  slot  type,  slant  slot  type,  and  turbulent  flow  type  with  stationary 
chamber).  Tests  were  conducted  on  a  homogeneous  flow  field,  quasi -homogeneous 
flow  field  (with  supporting  frame),  and  stability  state  circumferential 
inlet  distortion  flow  field.  The  experimental  work  was  conducted  at  a 
compressor  testing  platform  of  the  Shenyang  Aeroengine  Research  Institute. 

II.  STRUCTURE  CONFIGURATION  OF  CASING  TREATMENT 


Table  1  lists  four  types  of  structural  configuration  and  parameters  of 
casing  treatment.  The  physical  meaning  of  the  parameters  is  shown  in  Fig.  1. 

Table  1.  Structure  configurations  and  parameters  of  various  types  of 
casing  treatment. 
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Key:  (a)  Casing  types;  (b)  Axial  position;  (c)  Treatment  width; 

(d)  Included  angle  with  the  axial  line;  (e)  Ratio  of  opening  area; 
(f)  Radial  angle;  (g)  Height  of  stationary  chamber;  (h)  Others; 

(i)  Small  blade;  (j)  Circumferential  slot;  (k)  Slant  slot;  (1) 
Turbulent  flow  type  with  stationary  chamber;  (m)  Large  stationary 
chamber;  (n)  Slow  stationary  chamber;  (o)  Rear  axial  movement;  (p) 
Extruded  moving  blades;  (q)  Limited  to  core  of  moving  blades;  (r) 
Extruded  moving  blades;  (s)  Uniformly  distributed  255  blades  at 
circumference;  (t)  Uniformly  distributed  7  slots;  (u)  Uniformly 
distributed  360  blades  at  the  circumference;  (v)  Uniformly  distrib¬ 
uted  250  blades  at  the  circumference. 
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Fig.  1.  Treatment  casing  structural 
diagram  of  turbulent  flow  type  with 
stationary  chamber. 

Key:  (a)  Rotation;  (b)  Casing;  (c) 
Large  stationary  chamber;  (d)  Small 
stationary  chamber;  (e)  Axial  rear 
movement . 


III.  TESTING  RESUTLS  AND  ANALYSIS 


(A)  Comparison  of  increment  of  distortion  tolerance: 

The  increment  A  of  distortion  tolerance  is  calculated  according  to  the 
following  formula:  fRXW’  \  - 


In  the  equation,  R'  and  W'  are  the  converted  flows  and  pressure  ratio  of 
the  compressor  distortion  point  of  a  standard  casing. 

Figure  2  shows  the  compressor  characteristic  curves  of  a  standard  casing 
and  turbulent  flow  type  treatment  casing  with  stationary  chamber.  Generally, 
the  variation  rules  are  the  same  for  the  other  types  of  treatment  casing 
compressor  characteristic  curves.  Table  2  shows  the  increase  of  distortion 
tolerance  of  various  types  of  treatment  casing  compressors. 

It  is  apparent  from  Table  2  that  the  increment  of  compressor  distortion 
tolerance  of  quasi-uniform  flow  field  rated  revolution  speed  is  not  much. 


i 

i 


With  reduction  of  revolution  speed,  the  distortion  tolerance  increases 
appreciably;  the  greatest  increase  occurs  in  the  treatment  casing  with 
stationary  chamber. 


(a)  (b) 

Figure  2.  Compressor  performance  of  standard  casing  and  turbulent 
flow  type  treatment  casing  with  stationary  chamber:  (a)  Quasi -uniform 
ijilet  flow  field;  (b)  Circumferential  surge  inlet  flow  field  for 
D=0. 12. 

Key:  (1)  Kilograms  per  second;  (2)  Distortion  boundary  curves;  (3) 

Standard  distortion  boundary  curves;  (4)  Small  stationary  chamber; 

(5)  Standard  casing;  (6)  Large  stationary  chamber. 


Table  2.  Increments  of  compressor  distortion  tolerance  of  various  types 
of  casing  treatment  in  contrast  to  standard  casing. 


Key:  (a)  Inlet  flow  field;  (b)  Quasi^uniform  inlet  flow  field;  (c)  Circum¬ 
ferential  surge  inlet  flow  field  for  D«0. 12;  (d)  Casing  type;  (e)  Turbulent 
flow  type;  (f)  Large  stationary  chamber;  (g)  Small  stationary  chamber;  (h) 
Rearward  moving;  (i)  Small  blades;  (j)  Slant  slot;  (k)  Circumferential  slot. 
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Several  tests  of  inlet  surge  index  were  conducted  on  the  treatment 
casing  with  small  stationary  chamber.  At  high  revolution  speed,  the  inlet 
surge  index  is  increased,  and  the  increment  of  compressor  distortion 
tolerance  also  increases.  This  explains  that  casing  treatment  can  raise 
compressor  distortion  tolerance  at  high  revolution  speed. 

(B)  Comparison  of  the  maximum  efficiencies 

As  revealed  in  the  experimental  results,  for  various  types  of  inlet 
flow  field,  although  the  compressor  efficiency  curves  vary  for  various  types 
of  treatment  casing,  the  maximum  efficiency  values  approach  the  compressor 
efficiency  value  of  the  standard  casing. 

(C)  Comparison  of  compressor  outlet  pressure  field: 

In  a  circumferential  surge  flow  field  (D=0.12),  for  a  standard  casing 
at  the  rated  revolution  speed,  higher  non-uniformity  (D=8.7  percent)  of  circum¬ 
ferential  pressure  exists.  After  adopting  casing  treatment,  the  non -uniformity 
of  circumferential  pressures  at  compressor  outlet  was  significantly  reduced, 
being  more  pronounced  in  a  large  stationary  chamber  (D=6.7  percent).  However, 
these  only  improved  the  outlet  circumferential  pressure  field  of  the  treatment 
sector  without  variation  in  the  blade  root  region.  For  example,  for  the 
compressor  outlet  blade  tip  of  a  large  stationary  chamber,  D=8.7  percent  of 
the  standard  casing  was  reduced  to  4.9  percent. 

IV.  CONCLUSIONS 

(1)  Casing  treatment  can  effectively  increase  compressor  distortion 
tolerance  and  surge  margin;  the  results  are  related  to  structural  configurations. 
Size  of  the  stationary  chamber  and  axial  position  of  the  treatment  sector 

can  affect  compressor  performance.  This  explains  that  the  structural  configura¬ 
tions  and  parameters  should  be  rationally  selected. 

(2)  In  various  inlet  flow  fields,  the  maximum  values  of  different  casing 

treatment  compressors  approach  the  value  of  the  standard  casing. 
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(3)  In  a  surge  flow  field,  casing  treatment  can  improve  the  outlet  circum¬ 
ferential  pressure  field  of  the  blade  tip  region  at  high  revolution  speed. 

(4)  The  influences  on  increment  of  distortion  tolerance  of  several 
casing  treatments  are  the  same.  The  influence  is  higher  for  low  revolution 
speed  than  for  high  revolution  speed;  the  influence  of  circumferential 
surge  flow  field  is  greater  than  the  quasi-uniform  flow  field.  This  is 
the  natural  result  with  only  the  first  stage  casing  treatment  in  a  three 
stage  axial  compressor. 
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A  CALCULATION  METHOD  FOR  THE  ROCKET  ENGINE'S  OPTIMUM  THRUST  NOZZLE  CONTOUR 
DESIGN* 


Tian  Minhua 


Combining  the  characteristic  method  with  variation  principle  and  applying  it  to  su¬ 
personic  flow  field,  a  calculating  method  is  presented  for  the  optimum-thrust-nozzle 
contour  design.  This  method  can  be  used  to  design  nozzle  contours  with  various 
throat  curvature  radii  and  expansion  area  ratios.  It  is  especially  'suitable  for  the 
design  of  nozzle  contour  with  large  expansion  area  ratio. 

I.  CALCULATION  PRINCIPLE  OF  NOZZLE  CONTOUR  LINE  OF  OPTIMUM  THRUST 


Consider  an  ideal  two-dimensional  steady  equal  enthalpy  flow  without 
eddy.  The  flow  field  in  the  nozzles  is  axial  symmetric.  By  using  the 
cylindrical  surface  coordinate  system,  ABTE  in  Fig.  1  represents  the  inter¬ 
secting  curve  between  the  nozzle  contour  surface  and  the  meridian  surface. 
The  convergent  sector  AB  of  the  nozzle  throat  is  the  circular  arc  wall 
surface  with  Rj  as  the  radius  of  curvature.  The  initial  expansion  sector  BT 
at  the  lower  stream  of  the  throat  sector  is  the  circular  arc  wall  surface 
with  R2  as  the  radius  of  curvature.  The  author  and  his  colleagues  selected 
a  nozzle  contour  line  TE  to  have  the  maximum  thrust  produced  at  the  nozzle 
with  the  given  nozzle  length  L  and  environmental  pressure  pft.  This  contour 

*The  paper  was  read  at  the  Third  All  China  Thermophysics  Engineering  Confer¬ 
ence  at  Guilin  in  April  1980.  61 
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line  is  called  the  "optimal  thrust  nozzle  contour."  A  point  C  is  selected 
at  the  axial  line  of  the  nozzle.  Considering  a  control  surface  EC  passing 

A 

through  point  E  at  the  outlet  of  the  nozzle,  this  control  surface  EC  intersect 
at  point  D  with  the  right  characteristic  line  TH  passing  through  the  nozzle 
exit  point  E  (Fig.  1).  The  engine  thrust  R  can  be  derived  from  the  aerodynamic 
parameter  of  control  surface  EC. 

k  —  (‘  [  (r  —  ,.)  +  1  z-x'y  <0 

Jc  L  suj  <t>  J 

(p»  P,  v,  0  and  <P  indicate,  respectively,  the  pressure,  density,  velocity, 
gas  stream  angle,  and  angle  of  inclination  between  the  control  surface  and 
axis  of  the  gas  stream.) 


Fig.  1. 

Since  the  point  C  is  stationary,  the  nozzle  length  remains  constant. 
Therefore,  the  nozzle  length  after  C  also  remains  constant,  i.e., 

xt  —  *c  “  jc  ctg  4»Jy  —  constant  (2) 

In  addition,  because  of  continuity  of  flow,  the  flow  passing  through  the 
control  surface  should  be  equal  to  flow  G  passing  through  the  cross  sectional 
area  of  the  throat  sector;  i.e., 

ty 

Jc  cLn  4> 


G 


2  aydy  —  constant  (3) 


The  problem  is  then  summarized  into  the  situation  that  under  the 
condition  of  equations  (2)  and  (3),  the  thrust  R  of  the  nozzle  attains  the 
extreme  value;  i.e.,  the  arbitrary  function  equation  (1)  attains  the  extreme 
value  under  the  condition.  By  using  the  Lagrange  product  factor  [1],  the 
problem  of  extreme  value  under  the  condition  is  converted  into  the  unconditional 
extreme  value.  By  using  the  variable-difference  method,  the  thrust  R  can 
attain  the  maximum  value.  The  flow  parameters  on  control  surface  EC  should 
meet  the  various  following  conditions  [1]:  (Assume  that  the  external  pressure 
Pa=°) 


(1)  At  the  exit  terminal  point  E,  there  is 

rin2gB  -  (2/*Ml)  ctg  Of  (4) 

(2)  At  the  control  surface  ED,  there  is  the  left  characteristic  line: 

dy/dx  =  tg  (9+a) 

id  __  ctg  a  dp  +  ring «ing  _  Q  (5) 

dx  p  it  y  cos  (8  +  g) 

(3)  The  flow  parameters  on  ED  also  satisfy  the  two  following  relationship 
equations: 


Si* cat(0  —  g )/ co«g  “  Mi  eos(dt  —  gB)/ cosgB 
—  «  _  « 

-  M^j  ~  ria’gtgg  -■  ~ri&I0stggs 


In  the  equations. 

Si*  -  ( - i - V 

\  -  1  +  2/M1) 


(O 

(7) 


(4)  From  continuity  of  flow,  the  following  equation  is  apparently  true: 

Co "  5rfT7)  ’*>’*  -\\"  .)  <»> 

(k,  a,  and  M  indicate,  respectively,  the  adiabatic  index,  Mach  angle 
and  Mach  number  of  the  gas.) 


II.  GENERAL  DESCRIPTION  OF  CALCULATION  PROCESS 


In  nozzle  design,  generally  the  given  initial  conditions  are  and 
as  the  radii  of  the  throat  sector  of  the  nozzle,  gas  adiabatic  index  k,  and 
nozzle  expansion  area  ratio  Ag  in  order  to  obtain  the  optimal  thrust  nozzle 
contour.  For  nozzles  of  rocket  engine,  a  smaller  radius  of  curvature  at  the 
throat  sector  wall  surface  should  be  adopted  in  order  to  increase  thrust  and 
shorten  the  length.  The  author  and  his  colleagues  adopted  the  conformal 
curve  coordinate  method  in  manuscript  [2]  to  calculate  parameters  of  the 
transonic  initial  line  GK  (Fig.  1)  at  the  nozzle  throat  sector.  This  method 
is  simple  in  calculation  and  relatively  high  in  accuracy;  the  method  is 
adaptable  to  the  situation  of  a  relatively  small  radius  of  curvature  at  the 
throat  sector  wall  surface.  The  entire  supersonic  flow  field  calculation 
relied  on  a  characteristic  line  mesh  [3].  For  the  left  and  right  character¬ 
istic  line  equations,  difference  equations  of  two  stage  accuracy  are  adopted; 


i.e. , 


dy/ix  —  «g(0*±a*) 

A**ct*a*  ^  ±  Q 

v  >«*(£>•*<»•) 

(Parameters  with  indicate  the  average  values  for  Ax.) 


In  the  following,  the  derivation  method  (see  Fig.  1)  of  the  optimal 
thrust  nozzle  contour  TE  is  described.  From  the  characteristic  line  method, 
the  flow  parameters  from  the  wall  surface  to  the  right  characteristic  line 
T'H'  are  derived  item  by  item.  We  first  assume  a  from  equation  (4), 

9^  1  can  be  derived.  Then  a  point  D  is  sought  on  T'H'  to  satisfy  equation  (6). 
If  point  0  is  not  on  T’H’,  then  the  next  right  characteristic  line  is  derived. 

If  the  point  D  is  found,  y^1"1  can  be  determined  from  eq.  (7).  In  addition, 
parameters  of  various  points  on  ED  can  be  determined  from  eqs.  (5-7).  Moreover, 
by  using  eq.  (8),  we  can  determine  whether  T'  is  the  point  T  of  the  separation 
point.  If  eq.  (8)  is  not  true,  then  the  next  right  characteristic  line 
T'H'  and  point  D  are  derived  until  eq.  (8)  is  satisfied.  At  that  time,  T'H' 
is  the  TH  right  characteristic  line  required.  Again  based  on  y^1]  we  derive 
Ag1^  to  be  compared  with  the  given  Ag.  Usually,  Ag1^  is  not  equal  to  Ag.  By 
using  the  target  shooting  method,  a  suitable  Mg1^  can  be  selected.  We  repeated¬ 
ly  use  the  above  method  until  sufficiently  approaching  A_. 


Then,  the  right  characteristic  line  TD  and  parameters  on  control  surface 
ED  are  used  as  the  initial  condition,  based  on  TE  on  the  wall  surface  as  the 
streamline.  By  using  the  characteristic  line  method,  flow  field  parameters 
(Fig.  1)  of  the  TDE  triangular  region  can  be  calculated  to  obtain  the  optimal 
thrust  nozzle  contour  TE. 


III.  CALCULATION  RESULTS  AND  DISCUSSION 

In  the  following.  Table  1  lists  (for  different  area  expansion  ratio  Ag, 
and  throat  portion  expansion  sector  radius  R2)  parameters  of  separation  point 
(point  T)  and  the  outlet  terminal  point  (point  E)  parameters  for  the  optimal 
thrust  nozzle  wall  surface  contour.  In  the  table,  L  is  the  nozzle  length; 

CR  is  the  nozzle  thrust  coefficient. 


Table  1.  Rj  =  2,  k=1.23. 


^1 

*, 

*T 

Jr 

L 

*« 

Cm 

0.5 

0.277 

1.084 

33*36' 

8.399 

13*10' 

3.5147 

1.7605 

20 

1.0 

0.551 

1.165 

33*25' 

8.660 

13*08' 

3.5243 

1.7606 

1.5 

0.811 

1.238 

32*43' 

8.907 

13*06’ 

3.5326 

1.7609 

0.5 

R1 

1.096 

36*io* 

11*36' 

1.8261 

40 

1.0 

US 

1.189 

35*49* 

13.639 

11*34' 

4.0062 

1.8264 

1.5 

0.861 

1.272 

35*01' 

13.955 

11*33’ 

4.0128 

1.8265 

mm 

0.305 

1.104 

10*49” 

1.8586 

60 

mm 

0.604 

1.203 

IS 

10*48' 

4.2973 

1.8588 

1.5 

0.888 

1.291 

KSB 

Ha 

10*47' 

4.3031 

1.8589 

1  Considering  the  calculation  results,  for  a  certain  area  ratio,  the  smaller 

the  radius  R^  of  the  expansion  sector  of  the  nozzle  throat,  the  smaller  is 
,  the  nozzle  length  L.  In  order  to  shorten  the  nozzle  length,  we  can  appropriately 

reduce  R2;  however,  R2  cannot  be  extremely  small.  Because  the  smaller  the 
R2,  the  larger  is  the  expansion  angle  9^  of  the  initial  separation  point. 

When  R2  is  very  small,  the  gas  stream  at  the  throat  will  expand  exceedingly 
rapidly  to  increase  loss  of  the  gas  stream.  Therefore,  experiments  can  solve 
the  problem  of  what  is  the  optimal  value  for  a  small  R2  value.  Usually,  the 
radius  of  curvature  at  the  throat  wall  surface  can  be  taken  as  1$R^2  and 


O.S^R^l. 
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